Fundamental method and its hardware implementation for the generic prediction and analysis of multiple scattering of waves in particulate composites

ABSTRACT

This invention relates to a method and apparatus for the fundamental prediction, analysis, and parametric studies of the interaction of multiply scattered waves of any nature (acoustic, electromagnetic, elastic) with particulate composites (solid or fluid particulate phases included in a solid or fluid phase). It comprises of: a) a Prediction Engine  21  to predict the composite physical attributes from those of the constituent phases, b) an Inversion Engine  22  to estimate a subset of the attributes of the constituent phases and composite from the rest of those attributes, and c) an ergonomic Graphical Interface  23  to allow the user to easily set up complex simulation experiments to assist in parametric studies, analysis, and synthesis of real systems in which waves of any nature interact with composites. The Prediction Engine  22  is suitable for integration into existing acoustic or optical particle size analyzers to thereby extending substantially their concentration ranges.

[0001] Attached hereto there is a CD-R, hereby incorporated byreference, as an embodiment for the present invention. By followingindustry-standard installation procedures in a computer operating undera member of the MS-Windows® family of operating systems and a Pentium®II or higher equivalent processor, full disclosure for the currentinvention is enhanced, and its operability can be confirmed with actualexperimental data provided in the disk. A commercial product based onthis invention, to be called SCATTERER™ will be marketed in the nearfuture. The copyright owner has no objection to legal use andreproduction of this patent document and attached CD-R, as they appearin the Patent and Trademark Office records, reserving otherwise allrights.

BACKGROUND—FIELD OF INVENTION

[0002] This invention relates to a method and its hardwareimplementation for the fundamental prediction, analysis, and parametricstudies of the interaction of multiply scattered waves of any nature(acoustic, electromagnetic, elastic) with particulate composites (solidor fluid particulate phases included in a solid or fluid continuousphase). Any embodiment for the present invention constitutes a powerfulresearch and development tool in a vast variety of engineering andscientific fields dealing with wave propagation in composites. TheMETAMODEL™ Prediction Engine 21, a crucial component of this invention(FIG. 1), is suitable for integration into existing acoustic or opticalparticle size analyzers to thereby extending substantially theirconcentration ranges.

BACKGROUND—DISCUSSION OF PRIOR ART

[0003] Several systems to assist in the modeling, analysis, andparametric studies of physical systems are known in the prior art andavailable in the market. Some of them are general mathematical toolswith thousands of functions in calculus, linear algebra, ordinarydifferential equations, partial differential equations, Fourieranalysis, dynamic systems, etc. An incomplete list includes MATLAB™ (TheMathWorks, Inc), Maple V® (Waterloo Maple Software), SCIENTIFICWORKPLACE™ (TCI Software Research), MathCad™ (MathSoft, Inc.), andMathematica® (Wolfram Research, Inc.). These systems constitute, in ahierarchy of research and development tools, a next level up from e.g.the well-known IMSL FORTRAN/C library of scientific routines (VisualNumerics, Inc.). By becoming a little more specific in scope, they offera higher-level language to the user, reducing the need for proficiencyin computer programming and increasing the efficiency with which thescientist/engineer solves his/her problems. Some other systems are evenmore specific like FEMLAB® (COMSOL, Inc.) or MaX-1 (A visualElectromagnetics Platform, J. Wiley) or RealTime (A Two-DimensionalElectromagnetic Field Simulator, CRC Press). The last two systems areoverly specific simulators, as they cover only certain aspects ofelectromagnetic waves with no flexibility offered to modeling theirmultiple scattering in particulates. The first one, FEMLAB®, usingMATLAB™ as a software platform, offers to the user a less specifichigh-level language oriented towards the solution of partialdifferential equations using the finite element technique for problemsin physics and engineering. It allows the user to create a suitablemathematical model for the physical problem at hand from availablebuilding blocks (partial differential equations, functionalrelationships, etc.) in an ample variety of disciplines includingacoustics, fluid dynamics, structural mechanics, heat transfer, chemicalprocessing, and electromagnetics. After a short training, the user onlyneeds to be concerned with physics and engineering, and is required nomajor proficiency in computer programming and numerical intricacies.However, no pre-built modeling blocks for multiple scattering of wavesin particulate composites is available either.

[0004] All these systems in the prior-art are either exceedingly generalor overly specific to efficiently solve the multiply interleaved waveequations required to accurately predict the multiple scattering ofwaves in particulates. In order to take advantage of those existingtools, it would inevitably require an expert on both the physics and themathematics of the problem who—in addition—would have to command thenovel know-how to be disclosed in this specification, plus unacceptablecomputer time and human effort. A system tackling all those intricatephysico-mathematical aspects of the problem, and making them transparentto the regular physicist/engineer has been clearly in need for manyyears. The current invention provides such a tool.

[0005] Two distinct approaches have historically developed in trying tomathematically model scattering of waves with particulates. One ismicroscopic (fundamental) and called the analytical wave-scatteringtheory; the other is macroscopic (phenomenological) and called thetransport theory or radiative transfer theory ([22] A. Ishimaru, “WavePropagation and Scattering in Random Media”, Academic Press, 1978; [48]F. Alba et al, “Ultrasound Spectroscopy: A Sound Approach to Sizing ofConcentrated Particulates”, Handbook on Ultrasonic and DielectricCharacterization Techniques for Suspended Particulates, The AmericanCeramic Society, 1998). The former starts with the wave equation,rigorously predicts the scattering of a wave with a single particle inthe medium, and statistically describes the interaction between all thescattered fields around all the particles, to arrive to expected valuesfor the attributes of the composite in terms of the attributes of itsconstituents. It is a mathematically rigorous fundamental approach and,as a consequence, all the multiple scattering, diffraction, andinterference effects are contemplated with the broad validity andaccuracy range inherent to a first-principle technique. This is thephysico-mathematical background of the current invention. In the latter(macroscopic) approach, instead, the basic wave equation is not thestarting point. It deals directly with the transport of energy throughthe particulate composite, arriving to a basic differential equation(the equation of transfer). A basic assumption is that there is nocorrelation between the coexisting scattered fields and, therefore,their intensities can be added directly (as opposed to vectoriallyadding the fields). A crucial magnitude is the composite depth (particleconcentration by number times total cross-section of a single particletimes path length). The solution of the equation of transfer dependsupon this depth and the so-called phase function. The theory isheuristic as opposed to fundamental and, ergo, it lacks the mathematicalrigor and ample validity range of a first-principle approach. As itwould be expected, this approach is simpler and more tractable from thenumerical point of view. It has flourished during the last 100 years, ingreat part, due to the lack of sufficient cost-effective computationalpower to implement the fundamental approach. Nowadays, duplication ofcomputer power every year or so at lower costs, has made it possible toconsider full numerical implementation of the fundamental approach usingnon-expensive computer hardware.

[0006] Scientific literature on the analytical multiple-scatteringtheory is abundant but dispersed throughout a myriad of publications.Main authors who have significantly contributed in the last fifty yearsare L. L. Foldy, M. Lax, P. C. Waterman and R. Truell, V. Twersky, P.Lloyd and M. V. Berry, J. G. Fikioris and P. C. Waterman, A. K. Mal andS. K. Bose, N. C. Mathur and K. C. Yeh, V. N. Bringi V. V. Varadan, V.K. Varadan, Y. Ma, L. M. Schwartz, A. J. Devaney, L. Tsang and J. A.Kong, D. D. Phanord and N. E. Berger, de Daran, F. et al, etc. Work ofthese researchers and others are listed in this patent. The greatmajority of these publications are dedicated to specific problems inspecific fields of applications. A very few, however, recognizing thecommon underlying structural physico-mathematical core in all thesefields, have attempted to formulate generic equations for the predictionof either acoustical, electromagnetic, or elastic waves in particulatecomposites. The publication [37] “Multiple Scattering Theory forAcoustic, Electromagnetic, and Elastic Waves in Discrete Random Media”by V. V. Varadan, V. K. Varadan and Y. Ma (Multiple Scattering of Wavesin Random Media and Random Rough Surfaces, The Pennsylvania StateUniversity, 1985) belongs to this kind of general work and constitutesdirect theoretical prior-art to this invention as far as the PredictionEngine 21 is concerned (FIG. 1). Even though the approach described inthis prior-art and several others has a common theoretical backgroundwith the Prediction Engine 21, the latter constitutes a step forward ingenerality and accuracy, as the present specification will thoroughlydemonstrate.

[0007] There is finally a third approach in dealing with themultiple-scattering phenomenon that is more limited in scope: itconsists basically in attempting to suppress it by conceiving ad-hochardware designs and/or software schemes through which themultiple-scattered signal is either significantly attenuated in favor ofthe single-scattering component, or cross-correlation techniques areemployed to filter one component from the other. U.S. Pat. No.6,100,976, No. 5,956,139, and No. 4,380,817 belong to this group. Thisapproach is mentioned here for the sake of completeness even though,since no attempt to mathematically modeling multiple-scatteringphenomena is made, these patents do not truly constitute prior art tothe present invention.

[0008] Two important technical fields where the present invention meetsa long-term need are those of particle sizing using light scattering andacoustic spectroscopy. This is simply because, in order to indirectlymeasure a certain constituent attribute (e.g. particle sizedistribution) with a wave-based instrument, it is strictly necessary tobe able to predict—as a function of this constituent attribute—thatcomposite attribute that is directly being measured (e.g. scattering,attenuation, intensity). During the last 30 years or so, measuring ofparticle size in suspensions, emulsions, and aerosols using thescattering of light have developed to its mature state with a plethoraof commercial instruments available in the market. Nonetheless, inpractically all of those instruments a sine qua non condition foraccurate performance is that the particle concentration should bereduced before analysis to extremely low values (typically under 0.1%v/v)—when actual industrial concentrations are several orders ofmagnitude higher. One of the main reasons for such an operatingconstraint is that the mathematical model (e.g. Mie theory) internallyemployed by these light-scattering instruments is a single-scatteringmodel, i.e. a model that can predict accurately the composite attributesonly when the particles are well-separated between each other, so thatno interaction between contiguous electromagnetic fields occurs.Multiple scattering is simply not modeled.

[0009] A light-scattering instrument for aerosols (host medium isgaseous) that can work at higher concentrations (light transmissionsdown to 2% as opposed to typical values of about 60% to 80% in otherinstruments) is the one described by Harvill, T. L. and Holve, D. J. inU.S. Pat. No. 5,619,324, issued Apr. 8, 1997, in which correction formultiple scattering phenomena is conducted through a suitable adaptationof the theory described by E. D. Hirleman in [40] “Modeling of MultipleScattering Effects in Fraunhofer Diffraction Particle Size Analysis”,Particle Characterization 5, 57-65 (1988), and [41] “A General Solutionto the Inverse Near-Forward Scattering Particle Sizing Problem inMultiple Scattering Environments: Theory”, Proceedings of the 2^(nd)International Congress on Optical Particle Sizing, Mar. 5-8 1990, pp.159-168. With such a correction, this instrument extended theconcentration range by a factor of 5-10 (only up to about a few percentby volume). Hirleman's model was developed for the specific case ofsmall-angle scattering with the assumption that scattered fields addincoherently (on an intensity rather than amplitude basis). A “blackbox” approach is employed through defining a multiple-scatteringredistribution matrix that represents how optical energy incident in acertain direction is redistributed by the whole composite into anotherdirection. This matrix not only depends on the nature of medium andparticles with their size distribution and concentration, but also onthe geometry of the instrument as well as the optical depth of theparticulate system. In plain words, the matrix can only be determined(ad hoc) once the instrument (hardware) at hand is known. Hirleman usesalso a “successive orders of scattering” approach to express themultiple-scattering matrix in terms of the single-scattering one. Bothmatrices are global characterizations of the particulate system. Nofundamental modeling of the interaction of the wave with a singleparticle and interaction between scattered waves around particles isundertaken. By further assuming that a scattering event of order n canbe represented by the n-power of the single-scattering matrix applied tothe incident signal and summing all contributions, Hirleman arrives toan analytical relation between the single and multiple-scatteringmatrices. The black-box approach plus the dependence of the matrix upongeometry and extension of the medium gives the model a phenomenologicalcharacter as opposed to a fundamental (first principles) approach. Thenarrow applicability of a phenomenological approach is therefore ineffect. In U.S. Pat. No. 5,818,583, issued Oct. 6, 1998, E.Sevick-Muraca et al describe a method that applies the so-calleddiffusion model of light propagation to predict some aspects of multiplescattering of electromagnetic waves with suspensions (host medium is afluid). When the particle concentration is larger than about 1%, afurther approximation to the Radiative Transfer Theory can be madeleading to the Diffusion Approximation Theory ([22] A. Ishimaru, 1978).In this diffusion model, it is assumed that the diffusion intensity isscattered in all directions with an almost uniform angular distribution.In what the inventors describe as photon migration techniques, the lightintensity of the source is varied in time and the phase shift of thereceived signal for a number of different wavelengths is directlymeasured. The isotropic scattering coefficient is then calculated withthe diffusion model and, using this coefficient, the particle size andconcentration is determined inverting an integral equation that relatesthe particle size distribution with the scattering coefficient. Due tothe diffusion regime assumption, the model cannot be used at lowconcentration, because the presence of substantial multiple scatteringis a necessity for proper performance. Another limitation is that thereferred integral equation relating the scattering coefficient with thesize distribution assumes no interaction between particles occurs (thescattering coefficient is proportional to the particle concentration).This equation can only be valid when particles are far apart. As aresult, the model will only be accurate when multiple scattering isstrong enough for the diffusion approximation to hold but yet theparticles being far apart enough so no interaction occurs between themso equation 3 in column 9, line 55 of this prior-art patent can bevalid. These are opposite requirements in terms of particleconcentration. Consequent with this fact, the experimental data shown inthe patent correspond to concentrations around 1%—when actual industrialconcentrations are much higher than that. This is specifically discussedby the authors of the patent in [46] “Photon-Migration Measurement ofLatex Size Distribution in Concentrated Suspensions”, AIChE Journal,March 1997, Vol. 43, No. 3. In addition, due to the assumption of“almost uniform angular distribution of scattering”, validity of themodel will break down as the ratio between the size of the particles andthe wavelength increases. This is because the larger the particle sizeas compared to the wavelength, the stronger the directivity of thescattering in the forward direction is.

[0010] In U.S. Pat. No. 5,121,629, issued Jun. 16, 1992, I describe amethod and apparatus for measuring particle size distribution thatemploys a fundamental mathematical model for scattering of ultrasoundwaves in suspensions or emulsions (host medium is a fluid). The modelingaspects of this prior patent actually constitute the most directprior-art patent to the current invention as far as its application topredicting the attenuation of acoustic waves in suspensions isconcerned. Modeling details in this prior patent start in column 8, line29 up to column 12, line 46. As in the METAMODEL™ Prediction Engine 21of the current invention, modeling in this prior patent employs thefundamental wave scattering theory. It starts with the first-principledescription of the interaction of a wave with a single particle (theparticle/medium signature), and proceeds to statistically formulate thefield relations within an ensemble of particles of different sizes. Itends with an analytical equation that relates the propagation constantof the composite with the propagation constant of the host medium, thephysical properties of both medium and particles, as well as their sizedistribution and concentration. There are, however, three essentialdifferences with the current invention when applied to acoustics. Two ofthem are associated with the particle/medium signature for which thisprior-art patent uses the Epstein-Carhart-Allegra-Hawley (ECAH) theorywhile the present invention adds to this ECAH theory the following: a)modeling for possible viscoelastic behavior of both host medium andparticle and b) modeling for phenomena due to surface tension which isimportant in gaseous particles. The third paramount difference has to dowith the overlapping of contiguous fields: the prior patent arrives toequation 8 (column 10, line 60), which is the statistical relationobtained by Waterman and Truell in [9] “Multiple Scattering of Waves”,J. Math. Phys. Vol. 2, No.4, July-August, 1961, as a closed-formdispersion equation for the particulate composite. Such a simpleclosed-form equation was obtained under the assumption of statisticallyindependent point scatterers. In plain words, the expression “pointscatterers” refers to the abstract concept of a particle with nophysical dimensions but with the scattering properties of its actualsize. It was applied in order to simplify the mathematics and itsnumerical coding into computer software. As a result, model accuracydeteriorates as the size of the particles and/or frequency increase.Regarding the assumption of statistical independence between scattererspositions, it eventually breaks untrue as concentration increasesbecause of the inevitable spatial correlation between particles.Furthermore, P. Lloyd and M. V. Berry ([17] “Wave propagation through anassembly of spheres IV. Relations between different multiple scatteringtheories”, Proc. Phys. Soc., 1967, Vol. 91), proved there was an errorin the referred Waterman and Truell's formula “arising from anartificial stratification of the medium into thin slabs”. All thereferred drawbacks, though not damaging in many cases as demonstratedwith the good results disclosed in this prior patent, would insteaddefinitely contribute to destroy the accuracy of a generic model (validfor waves of any nature) if the broadest range of materials, particlesize, concentration, and wave frequencies is aimed at. In fact, it willbe demonstrated when the Prediction Engine 21 is described in detailthat, if scatterers are treated as what they are: finite-size inclusionsin the host medium, and the spatial correlation between them is notneglected, it is not possible any longer to arrive at a simpleclosed-form expression for the dispersion equation. Mathematics andcomputer programming complexities increase considerably—but predictionaccuracy and range of applicability improve drastically.

[0011] In U.S. Pat. No. 6,109,098, issued Aug. 29, 2000, A. Dukhin andP. Goetz describe a device that combines acoustic and electroacousticspectrometry to measure particle size distribution and zeta potentialfor concentrated dispersed systems. Focusing only on the mathematicalmodeling part of this prior patent, they use a macroscopicphenomenological approach to modeling the interaction between the soundwave and the concentrated suspension (host medium is a fluid).Specifically, they use the so-called “coupled phase”” and “cell” modelsapplying the long-wave regime condition, i.e. the particle size has tobe much smaller than the sound wavelength. As a result, their inventiondelivers useful results only for particles with sizes smaller than 10μm. It is worth noting that the authors categorically reject beforehandthe very fundamental approach adopted in the present invention whendiscussing U.S. Pat. No. 5,121,629 by saying in column 4, line 55: “. .. He (Alba) suggests using EKAH (they mean ECAH) theory in combinationwith a multiple scattering approach This is not adequate because eventhe multiple scattering approach requires a single particle acousticfield which is known only for a single particle in infinite media . . .” FIG. 3 through FIG. 11 shown in the current specification willirrefutably demonstrate that the Prediction Engine 21 of the presentinvention has accomplished what U.S. Pat. No. 6,109,098 says it wasunfeasible.

[0012] In U.S. Pat. No. 6,119,510, issued Sep. 19, 2000, M. Carasso. etal describe an “improved process for determining the characteristics ofdispersed particles”. The authors contemplate the use of their processwith acoustic as well as light waves. Focusing again only in themathematical modeling, they describe a “computationally improvedAllegra-Hawley Model” for acoustics and the original Mie Model for lightwaves. As already stated in depth, those models are single-scatteringmodels that will inevitably fail as the concentration increases(typically around 10% for acoustics and 1% for optics). No attemptwhatsoever is made to model multiple-scattering phenomena—not even byemploying equation 8 (column 10, line 60) of U.S. Pat. No 5,121,629 thatis the closest established prior-art per their own admission. Theauthors simply say: “Although this predicted spectrum is anapproximation (owing to non-linearities induced by multiple-scatteringphenomena) it provides desirable results in this invention”. The meaningof the word “desirable” is nonetheless never tested because the acousticexperimental data shown in this prior patent correspond toconcentrations of 1% and 5% for which multiple scattering is negligible.Similarly, the optical turbidity data correspond to a concentration ofonly 1%.

[0013] In summary, there are no research and development tools availablefor the generic prediction, analysis, and parametric studies of theinteraction between generic waves (acoustic, electromagnetic, elastic)and particulate composites. The current invention meets such a long-termneed by providing a powerful research and development tool that makesall theoretical physics and numerical intricacies transparent to thehuman operator. A fundamental mathematical model is a key component ofsuch a tool. The lack of a complete fundamental mathematical model forthe prediction of multiple scattering of waves in particulates has, inaddition, limited the concentration range of currently existingwave-based instruments. All such instruments which operation is based onthe interaction of any type of wave traveling through suspensions,emulsions or solids containing inclusions, could benefit considerablyfrom the prediction engine of this invention. Because the METAMODEL™Prediction Engine 21 is a separate module connected to the other modulesof this invention through industry-standard interfaces, it can bedirectly integrated to those—already existing—instruments extendingtheir concentration ranges several orders of magnitude. Further objectsand advantages of my invention will become apparent from a considerationof the ensuing description.

SUMMARY OF THE INVENTION

[0014] This invention comprises of three crucial components (FIG. 1): a)a fundamental mathematical model, with its corresponding numericalalgorithms to predict the composite physical attributes from thephysical attributes of the constituent phases. I call this component theMETAMODEL™ Prediction Engine 21, b) a set of optimization algorithms tonumerically invert the aforesaid mathematical model so as to estimate asubset of the attributes of the constituent phases and composite fromthe knowledge of the rest of the attributes of the constituent phasesand composite. I call this component the METAMODEL™ Inversion Engine 22;and c) an ergonomic graphical interface between the prediction andinversion engines and the final user so he/she can easily set up complexsimulation experiments to assist in parametric studies, analysis, andsynthesis of real systems in which waves of any nature interact withcomposites. I call this component the SCATTERER™ Graphical Interface 23.The Prediction Engine 21 uses a fundamental approach that considers thecomposite as a random medium and consequently describes it statisticallyas indicated in block 24 of FIG. 2. The particle/host medium signatureis deterministically modeled from the knowledge of the physicalproperties of particle and host medium and the physics associated withthe nature of the problem (acoustical, electromagnetic, elastic) asexpressed in block 25. Particle spatial and attribute correlations aremodeled as indicated in block 26. Stochastic field equations, assuming aparticle is fixed in space with its physical attributes known, areestablished as represented in block 27. Further averaging over thephysical attributes of the fixed particle is performed as expressed inblock 28, arriving to the dispersion equations that are solved todetermine the propagation constant of the random medium as indicated inblock 29. All physical attributes of the composite are thus determinedfrom the composite propagation constant and the physical attributes ofthe constituent phases as depicted in block 30. The graphical interface23 is implemented following de facto industry standards for graphicaluser interfaces. Hierarchies of windows, icons, options and commandbuttons, checkbox controls, context menus, image controls, tool-tips,on-line help and the like are used to produce an ergonomic efficientinteraction between this invention and the human operator. Ageneral-purpose computer 19 and data-carrier media 20 are needed for thefull operability of the invention.

DESCRIPTION OF THE DRAWINGS

[0015]FIG. 1 schematically depicts main components of the invention withtheir structural relationship.

[0016]FIG. 2 schematically displays main modeling steps in theMETAMODEL™ prediction Engine.

[0017]FIG. 3 validates the prediction engine for an extremely finedispersion of silica particles (19 nm) used in the semiconductorindustry for polishing, at a particle concentration of C_(v)=11.4%

[0018]FIG. 4 validates the prediction engine for an extremely finedispersion of silica particles (19 nm) used in the semiconductorindustry for polishing at C_(v)=17.1%

[0019]FIG. 5 validates the prediction engine for a dispersion ofstandard silica particles (100 nm±30 nm) at C_(v)=10.0%

[0020]FIG. 6 validates the prediction engine for a dispersion ofstandard silica particles (100 nm±30 nm) at C_(v)=30.0%

[0021]FIG. 7 validates the prediction engine for a dispersion ofstandard silica particles (300 nm±30 nm) at C_(v)=16.0%

[0022]FIG. 8 validates the prediction engine for a dispersion ofstandard silica particles (300 nm±30 nm) at C_(v)=24.0%

[0023]FIG. 9 validates the prediction engine for a dispersion ofstandard silica particles (450 nm±30 nm) at C_(v)=10.0%

[0024]FIG. 10 validates the prediction engine for a dispersion ofstandard silica particles (450 nm±30 nm) at C_(v)=30.0%

[0025]FIG. 11 validates the prediction engine for a broad distributionof soy oil droplets in water (0.06 μm to 2 μn) at C_(v)=10.6%

DETAILED DESCRIPTION/OPERATION OF THE INVENTION

[0026] Detailed Description/Operation of the METAMODEL™ Engines 21 and22

[0027] An essential part of this invention is a generic fundamentalmathematical model to predict the interaction between a generic wave(acoustic, electromagnetic, elastic) and a dense generic particulatesystem (fluid or solid particles included in fluid or solid media). Twonumerical engines operate around this mathematical model. The PredictionEngine 21 starts with the operating attribute vector o, the attributevector m of the host medium, and the attribute vector a of the particlesto predict the attribute vector c of the composite (FIG. 1). Vector ocontains operating parameters like frequency, temperature, scatteringangle, etc. Vector m contains all the physical properties for thesuspending medium. Vector a comprises of all the physical properties forthe material of the particles plus the particle size distribution andconcentration. Some of these properties can be frequency-dependent (e.g.for non-Newtonian fluid or viscoelastic solids, or the refractive indexin the optical case). Properties may also be temperature-dependent. Thecomponents of vector c are the physical attributes of the composite(attenuation, velocity, intensity, phase shift, elastic constants,etc.). These attributes are in general functions of frequency,temperature, scattering angle, etc. As a clarifying example, for theacoustic case treated in U.S. Pat. No. 5,121,629 the knowledge of theseven disclosed physical properties for both the suspending medium andsuspended particles plus the particle size distribution and particleconcentration, made it possible to uniquely predict the attenuationcoefficient of the ultrasound wave as a function of frequency. TheInversion Engine 22 handles the inverse problem, i.e. knowing some ofthe attributes of the composite and constituents, it estimates thevalues of the missing attributes of the constituents and composite (FIG.1). The symbol P{.} stands for a subset of all properties (the partialknowledge set). The symbol {overscore (P)}{.} stands for the missingsubset (set complement). As a concrete example, most of the referredpatents describe particle size analyzers, viz. instruments that estimatethe particle size distribution (a constituent attribute) from themeasurement of a composite attribute (attenuation spectrum for U.S. Pat.No. 5,121,629 and U.S. Pat. No. 6,109,098, light intensity as a functionof angle for U.S. Pat. No. 5,619,324, and phase shift as a function ofwavelength for U.S. Pat. No. 5,818,583) and the knowledge of the rest ofthe attributes. All the inherent numerical intricacies associated withthe inversion of integral equations are well known by those skilled inthe art. Non-linear optimization algorithms using Powell or relatedtechniques in multi-dimensional space (W. H. Press et al, “NumericalRecipes, The Art of Scientific Computing”, Cambridge Press, 1988) havebeen successfully employed for the specific acoustic problem describedin U.S. Pat. No. 5,121,629. The present invention, using similar butmore powerful and elaborate numerical techniques, provides for ageneralized direct/invert algorithmic scheme allowing the human operatorto manipulate the mathematical model in any desired direction specifyingthose parameters that are known as well as those that are not known andare to be estimated. As an example, the engineer/scientist could set upthe current invention to use measured attenuation spectral data,particle size, and concentration as inputs to estimate the elasticconstants of the particles.

[0028] In the attached embodiment for this invention, both numericalengines 21 and 22 are programmed in Compaq Visual FORTRAN as aDynamically Linked Library (DLL) called by the Graphical Interface 23that is programmed in Microsoft Visual Basic® for Windows. Any otherlanguage or combination of languages interacting with a different likingtechnology, e.g. Component Object Model (COM) server technology, couldobviously be used to embody the current invention.

[0029] In this generic modeling problem I start considering a randomspatial distribution of N→∞ particles (scatterers) with arbitraryphysical attributes (thermodynamical, size, shape, orientation,acoustic, elastic, electromagnetic, transport, etc.) suspended orincluded in an otherwise continuous medium with an arbitraryconcentration. The particular position and attributes of each scattererare unknown but the ensemble of them can be statistically describedthrough a probability density function p(r ₁, . . . ,r _(i), . . . ,r_(N);a ₁, . . . ,a _(i), . . . ,a _(N)). The symbols r _(i) and a _(i)represent the position and attribute vectors for particle irespectively.

[0030] In plain words: p(r ₁, . . . ,r _(i), . . . ,r _(N);a ₁, . . . ,a_(i), . . . ,a _(N))dv₁ . . . dv_(i) . . . dv_(N)da ₁ . . . da _(i) . .. da _(N) is the joint probability of finding one scatterer in thevolume dv₁ (centered at r ₁) with attribute vector in the neighborhoodof a ₁ (|a−a ₁|≦ε), another scatterer in dv₂ with attribute vector inthe neighborhood of a ₂, etc. The symbol dv_(i) stands for the volumedifferential element in R³ with R being the set of real numbers. CallingnA the number of attributes for each particle, note that the symbol da_(i) stands not for the differential of the vector a _(i), but for thehyper-volume differential element in R^(nA). One of the components of a_(i) is the size of particle i that will be denoted as d_(i). Thesimplest practical case commonly considered in the priormultiple-scattering literature is that in which attributes of allscatterers are identical. A few authors considered as well the case inwhich there is a distribution of sizes throughout the particlepopulation.

[0031] A generic wave incites the above random medium. As it is wellknown from the prior art, assuming the particle concentration is lowenough, the local exciting wave at each scattering center can beconsidered to be simply the global incident wave (Born's approximation).This is the well-known single-scattering regime where the behavior ofthe ensemble is simply proportional to the behavior of the individualparticles and the number of them per unit volume. For the acousticalcase, this incident wave is normally a longitudinal wave(compressional), i.e., its displacement vector is parallel to thedirection of propagation of the wave. At each scattering center thislongitudinal wave converts into multiple modes inside as well as outsidethe particle: two compressional waves (refracted and reflected), twothermal longitudinal waves (inside and outside), and twoviscous-inertial transversal waves (inside and outside). At lowconcentrations these new modes extinguish in the close neighborhood ofeach particle before interacting with other particles. As theconcentration increases, these local waves around the particles overlapand multiple-scattering effects become significant. The detected wave atthe receiver is the remaining propagating compressional mode. For theelectromagnetic case, the incident wave is transversal, i.e., the planeof vibration of the electrical and magnetic fields is perpendicular tothe direction of propagation of the wave. There are normally noconversion modes at each scatterer. The elastodynamic case is the mostcomplex of all, as its vectorial nature is present in both the incidentwave as well as in the scattering process that occurs at each particle.The incident wave can have both longitudinal and transversal componentsand each one of them converts—at each scattering center—intolongitudinal and transversal modes.

[0032] The incident field is assumed to be that of a monochromatic wavewithout loss of generality. The factor e^(-jω) is consequently implicitin all equations. Except for the incident field, all fields in theparticulate system depend upon the spatial/attribute configuration ofscatterers {r ₁, . . . ,r _(i), . . . ,r _(N);a ₁, . . . ,a _(i), . . .,a _(N)}. This lengthy notation will be maintained throughout thepresent description for the sake of an enabling disclosure. I call theincident field at observation point r by I(r), and the field thatexcites particle i (located at point r _(i) with attribute vector a_(i)) at observation point r by E(r/r _(i);a _(i):{r ₁, . . . ,r _(i), .. . ,r _(N);a ₁, . . . ,a _(i), . . . ,a _(N)}). Similarly, calling thefield scattered by particle j (located at r _(j) with attribute a _(j))when particle i with attribute a _(i) is located at r _(i) by NS(r/r_(j);r _(i);a _(j),a _(i):{r ₁, . . . ,r _(i), . . . ,r _(j), . . . ,r_(N); a ₁, . . . ,a _(i), . . . ,a _(j), . . . ,a _(N)}), the followingimportant generic equation is valid: $\begin{matrix}{{\underset{\_}{E}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{{\underset{\_}{a}}_{i}:\left\{ {{{\underset{\_}{r}}_{1}{,\quad}\quad \ldots {\quad,\quad}\quad {\underset{\_}{r}}_{i},\quad \ldots {\quad,}\quad {\underset{\_}{\quad r}}_{N}};{{\underset{\_}{a}}_{1},\quad \ldots \quad {,\quad}\quad {\underset{\_}{a}}_{i}{,\quad}\quad \ldots \quad {,\quad}\quad {\underset{\_}{\quad a}}_{N}}} \right\}}} \right)} = {{\underset{\_}{I}\left( \underset{\_}{r} \right)} + {\sum\limits_{j = {{1,\quad j} \neq i}}^{N}\quad {\underset{\_}{S}\left( {{{\underset{\_}{r}/{\underset{\_}{r}}_{j}}{\underset{\_}{r}}_{i}};{{\underset{\_}{a}}_{j}{{\underset{\_}{a}}_{i}:\left\{ {{{\underset{\_}{r}}_{1}{,\quad}\quad \ldots {\quad,}\quad {\underset{\_}{r}}_{i}{,\quad}\quad \ldots \quad {,\quad}{\underset{\_}{r}}_{j},\quad \ldots \quad,\quad {\underset{\_}{r}}_{N}};{{\underset{\_}{a}}_{1}{,\quad}\quad \ldots \quad {,\quad}\quad {\underset{\_}{a}}_{i}{,\quad}\quad \ldots \quad {,\quad}\quad {\underset{\_}{a}}_{j}{,\quad}\quad \ldots \quad {,\quad}\quad {\underset{\_}{a}}_{N}}} \right\}}}} \right)}}}} & (1)\end{matrix}$

[0033] Following Foldy's [2] approach to handle the statistical natureof the physical system under consideration, I define the expected(average) value of a field at the observation point r as theconfigurational average over all possible spatial/attributedistributions of the scatterers: $\begin{matrix}{{{\langle{\underset{\_}{F}\left( \underset{\_}{r} \right)}\rangle} = {\int_{V}^{\quad}{\ldots {\int_{V}^{\quad}{\int_{D}^{\quad}{\ldots {\int_{D}^{\quad}{{\underset{\_}{F}\left( \quad {\underset{\_}{r}:\left\{ {{{\underset{\_}{r}}_{1}{\quad,\quad}\quad \ldots \quad,\quad {\underset{\_}{r}}_{i}{\quad,}\quad \ldots {\quad,\quad}\quad {\underset{\_}{r}}_{N}};{{\underset{\_}{a}}_{1}{\quad,}\quad \ldots \quad,\quad {\underset{\_}{a}}_{i}{\quad,}\quad \ldots \quad {\underset{\_}{a}}_{N}}} \right\}} \right)}{p\left( {{\underset{\_}{r}}_{1}{\quad,}\quad \ldots \quad,\quad {\underset{\_}{r}}_{i}\quad,\quad \ldots \quad {{\underset{\_}{r}}_{N};{{\underset{\_}{a}}_{1}{\quad,}\quad \ldots {\quad,}\quad {\underset{\_}{a}}_{i\quad}{,\quad}\quad \ldots {\quad,}\quad \ldots {\quad,}\quad {\underset{\_}{a}}_{N}}}} \right)}{v_{1\quad}}\ldots \quad {v_{i}}\ldots \quad {v_{N}}{\quad {\underset{\_}{a}}_{1\quad}}\ldots \quad {{\underset{\_}{a}}_{\quad_{i}\quad}}\ldots \quad {\underset{\_}{a}}_{N}\quad {\forall{\underset{\_}{F} \in \left( {\underset{\_}{I},\quad \underset{\_}{E},\quad \underset{\_}{S}} \right\}}}}}}}}}}}\quad} & (2)\end{matrix}$

[0034] The domain of integration V is the volume accessible to thescatterers, i.e., a subset of R³. The domain D of the attributes is asubset of R^(nA). Because of the random nature of the system all fieldscan be thus expressed as:

F ( r )=< F ( r )>+ F ₁ ( r )  (3)

[0035] Symbols <F(r)> and F ₁ (r) represent the coherent and incoherentfields respectively. The incoherent field is the truly random part ofthe total field. Through the introduction of conditional probabilities,the expected value of a field observed at r when a particle is assumedto be at r _(i) with attribute vector a _(i) (first order conditionalaverage) is: $\begin{matrix}{{\langle{\underset{\_}{F}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\rangle} = {\int_{Vi}^{\quad}{\ldots^{\prime}\ldots {\int_{Vi}^{\quad}{\int_{D}^{\quad}{\ldots^{\prime}\ldots {\int_{D}^{\quad}{{\underset{\_}{F}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{{\underset{\_}{a}}_{i}:\left\{ {{{\underset{\_}{r}}_{1}{,\quad}\quad \ldots {\quad,\quad}\quad {\underset{\_}{r}}_{i}{\quad,}\quad \ldots \quad,\quad {\underset{\_}{r}}_{N}};{{\underset{\_}{a}}_{1}\quad {,\quad}\quad \ldots \quad,\quad {\underset{\_}{a}}_{i\quad},\quad \ldots \quad,\quad {\underset{\_}{a}}_{N}}} \right\}}} \right)}{p\left( {{{\underset{\_}{r}}_{1}{\quad,}\quad \ldots \quad {,\quad}^{\prime}\ldots \quad {\underset{\_}{r}}_{N}};{{\underset{\_}{a}}_{1\quad},\quad \ldots \quad {,\quad}^{\prime},\quad \ldots {\quad,}\quad {{\underset{\_}{a}}_{N}/{\underset{\_}{r}}_{i}}};{\underset{\_}{a}}_{i}} \right)}{{v\quad}_{1\quad}},\quad \ldots^{\prime}\ldots \quad {v_{N}}\quad {\quad {\underset{\_}{a}}_{1}}\quad \ldots^{\prime}\ldots \quad {{\underset{\_}{a}}_{N}}\quad {\forall{\underset{\_}{F} \in \left\{ {\underset{\_}{I},\quad \underset{\_}{E},\quad \underset{\_}{S}} \right\}}}}}}}}}}} & (4)\end{matrix}$

[0036] The primes in the previous formula mean that r _(i) and a _(i)are not varied throughout the averaging process as clearly expressed bythe conditional probability. The domain of integration Vi is V minus asphere of diameter d_(i). For the case of assuming particles i and jrespectively at locations r _(i) and r _(j) with properties a _(i) and a_(j) (second order conditional average) it is obtained: $\begin{matrix}{\left. {{\langle{\underset{\_}{F}\left( {{{\underset{\_}{r}/{\underset{\_}{r}}_{i}}{,\quad}\quad {\underset{\_}{r}}_{j}};{{\underset{\_}{a}}_{i}{,\quad}\quad {\underset{\_}{\quad a}}_{j}}} \right)}\rangle} = {{\int_{Vij}^{\quad}{\ldots^{\prime}\ldots^{\prime}\ldots \quad {\int_{Vij}^{\quad}{\int_{D}^{\quad}{\ldots^{\prime}\ldots^{\prime}\ldots {\int_{D}^{\quad}{\underset{\_}{F}\left( {{{\underset{\_}{r}/{\underset{\_}{r}}_{i}},\quad {\underset{\_}{r}}_{j}};{{\underset{\_}{a}}_{i},\quad {\underset{\_}{a}}_{j}}} \right)}}}}}}}:\left\{ {{{\underset{\_}{r}}_{1},\quad \ldots \quad,\quad {\underset{\_}{r}}_{i},\quad \ldots {\quad,}\quad {\underset{\_}{r}}_{j},\quad \ldots \quad,\quad {\underset{\_}{r}}_{N}};{{\underset{\_}{a}}_{1},\quad \ldots \quad,\quad {\underset{\_}{a}}_{i},\quad \ldots {\quad,}\quad {\underset{\_}{a}}_{j},\quad \ldots \quad,\quad {\underset{\_}{a}}_{N}}} \right\}}} \right){p\left( {{{\underset{\_}{r}}_{1}\ldots {,\quad}^{\prime},\quad \ldots {,\quad}^{\prime},\ldots \quad {\underset{\_}{r}}_{N}};{{\underset{\_}{a}}_{1},\ldots ,^{\prime},\ldots ,^{\prime},\ldots {,\quad}{{\underset{\_}{a}}_{N}/{\underset{\_}{r}}_{i}},\quad {\underset{\_}{r}}_{j}};{{\underset{\_}{a}}_{i},\quad {\underset{\_}{a}}_{j}}} \right)}\quad {v_{1}}\ldots^{\prime}\ldots^{\prime}\ldots \quad {v_{N}}\quad {{\underset{\_}{a}}_{1}}\ldots^{\prime}\ldots^{\prime}\ldots \quad {{\underset{\_}{a}}_{N}}\quad {\forall{\underset{\_}{F} \in \left\{ {\underset{\_}{I},\quad \underset{\_}{E},\quad \underset{\_}{S}} \right\}}}} & (5)\end{matrix}$

[0037] Without further explanation, statistics of higher order arestraightforwardly defined. The global configurational average <F(r)> canthen be expressed in terms of the first partial average as follows:$\begin{matrix}{{\langle{\underset{\_}{F}\left( \underset{\_}{r} \right)}\rangle} = {\int_{V}^{\quad}{\int_{D}^{\quad}{{\langle{\underset{\_}{F}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\rangle}{p\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\quad {{\underset{\_}{a}}_{i}}\quad {v_{i}}\quad {\forall{\underset{\_}{F} \in \left\{ {\underset{\_}{I},\quad \underset{\_}{E},\quad \underset{\_}{S}} \right\}}}}}}} & (6)\end{matrix}$

[0038] Similarly, the first partial average can be expressed in terms ofthe second as: $\begin{matrix}{{\langle{\underset{\_}{F}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\rangle} = {{\int_{Vi}^{\quad}{\int_{D}^{\quad}{{\langle{\underset{\_}{F}\left( {{{\underset{\_}{r}/{\underset{\_}{r}}_{i}},\quad {\underset{\_}{r}}_{j}};{{\underset{\_}{a}}_{i},\quad {\underset{\_}{a}}_{j}}} \right)}\rangle}{p\left( {{\underset{\_}{r}}_{j};{{\underset{\_}{a}}_{j}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\quad {{\underset{\_}{a}}_{j}}\quad {v_{j}}}}} = {\int_{Vi}^{\quad}{\int_{D}^{\quad}{{\langle{\underset{\_}{F}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{j}};{\underset{\_}{a}}_{j}} \right)}\rangle}{p\left( {{\underset{\_}{r}}_{j};{{\underset{\_}{a}}_{j}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\quad {{\underset{\_}{a}}_{j}}\quad {v_{j}}\quad {\forall{\underset{\_}{F} \in \left\{ {\underset{\_}{I},\quad \underset{\_}{E},\quad \underset{\_}{S}} \right\}}}}}}}} & (7)\end{matrix}$

[0039] In general then, the partial average of order k can be expressedin terms of the partial average of order k+1 and the correspondingconditional probability. Pursuing this recurrence throughout the wholepopulation is obviously impractical. In order to break the recurrence, Iused in the above second equality a generalization of the so-called“quasi-crystalline” approximation introduced first by Lax [3,4]. Itconsists in assuming that <F(r/r _(i),r _(j);a _(i),a _(j)(>=<F(r/r_(j);a _(j))>. Breaking the recurrence at a higher statistical level isalso contemplated in this invention but it would complicate considerablymore the mathematical modeling process as well as its disclosure in thisspecification with no significant practical advantage.

[0040] Taking now the configurational average of both terms in equation1 and using equation 7, it is thus obtained: $\begin{matrix}{{\langle{\underset{\_}{E}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\rangle} = {{\underset{\_}{I}\left( \underset{\_}{r} \right)} + {\sum\limits_{j = {{1,\quad j} \neq 1}}^{N}\quad {\int_{Vi}^{\quad}{\int_{D}^{\quad}{{\langle{\underset{\_}{S}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{j}};{\underset{\_}{a}}_{j}} \right)}\rangle}{p\left( {{\underset{\_}{r}}_{j};{{\underset{\_}{a}}_{j}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\quad {{\underset{\_}{a}}_{j}}\quad {v_{j}}}}}}}} & (8)\end{matrix}$

[0041] Let us now define n(r _(j),a _(j))da _(j)dv_(j) as the number ofscatterers in the neighborhood of point r _(j) and with attribute vectora in the neighborhood of a _(j). Let us further define p(r _(j);a_(j))da _(j)dv_(j)=n(r _(j),a _(j))da _(j)dv_(j)/N, i.e., the ratiobetween the number of particles with attribute in the neighborhood of a_(k) and location in the neighborhood of r _(j), and the total number ofparticles in the ensemble. Assuming there is no global correlationbetween position and attribute, it is possible to simplify the jointprobability as p(r _(j);a _(j))=n(r _(j))n(a _(j))/N²=n(r _(j))f_(n)(a_(j))/N. We have introduced the function f_(n)(a _(j)) as the densitydistribution by number for the attributes of particle j in theneighborhood of r _(j). When the attribute is particle size, f_(n)(a_(j)) is the standard particle size distribution by number but local tothe neighborhood of r _(j). The non-segregation assumption also impliesf_(n) (a _(j))=f_(n) (a _(i))=f_(n) (a) which is now the standardensemble particle size distribution. Along similar intuitive lines, theconditional probability can be defined as:

p( r _(j) ;a _(j) /r _(i) ;a _(i))={n( r _(j) ,a _(j))/(N−1)}g( r _(j),a _(j) ,r _(i) ,a _(i))   (9)

[0042] The function g(r _(j),a _(j),r _(i),a _(i)) is a measure of thelocal correlation between the positions and attributes of the particles.It basically depends on the nature and range of the interparticleforces. A very successful approximation commonly used in statisticalmechanics (the so-called Pair Correlation Function) consists in assumingparticles are impenetrable and that the correlation is only a functionof their distance r=|r _(j)−r _(i)| and the sum of their radiiD=(d_(j)+d_(i))/2:

g( r _(j) ,a _(j) ,r _(i) ,a _(i))=g(r,D)=H(r−D)CF(r)   (10)

[0043] The function H(x) is the Heaviside step function which value is 0for x<0 and 1 for x≧0. When CF(r)=1 it is the simplest approximation tocorrelation: as long as particles do not interpenetrate, given thatparticle i is located at r _(i), the chance of particle j to be locatedat r _(j) with attribute a _(j) is governed exclusively by p(r _(j),a_(j)). This is called the “hole correction” in the prior art [13]. Ingeneral, the function CF(r) is oscillatory approaching unity as thedistance between the two particles increases simply because as twoparticles separate away, their spatial correlation disappear. There isan extended literature regarding the calculation of CF(r). The mostsuccessful and well-known form is the so-called Percus-Yevickapproximation from the theory of classical fluids. Other approaches arethe so-called hypernetted chain equation, the Virial expansion, theself-consistent approximation, etc. A powerful technique recentlydeveloped is based on Monte-Carlo simulation of a large number ofscatterer spatial configurations. The current invention, in itsgraphical user interface, provides for most of these well-known formsfor CF(r) as set-up options. Another powerful option provided by thisinvention is to allow for this oscillatory function to be characterizedby three parameters (amplitude, frequency, and decay) and let the userexperiment with them in trying to accurately modeling the unique spatialcorrelation characteristics of the problem at hand. The followingrelation is then valid: $\begin{matrix}{{\langle{\underset{\_}{E}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\rangle} = {{\underset{\_}{I}\left( \underset{\_}{r} \right)} + {\int_{Vi}^{\quad}{\int_{D}^{\quad}{{\langle{\underset{\_}{S}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{j}};{\underset{\_}{a}}_{j}} \right)}\rangle}{n\left( {{\underset{\_}{r}}_{j},\quad {\underset{\_}{a}}_{j}} \right)}{g\left( {{\underset{\_}{r}}_{j},\quad {\underset{\_}{a}}_{j},\quad {\underset{\_}{r}}_{i},\quad {\underset{\_}{a}}_{i}} \right)}\quad {{\underset{\_}{a}}_{j}}\quad {v_{j}}}}}}} & (11)\end{matrix}$

[0044] The field scattered by a single particle has been separatelystudied in the prior literature for all acoustic, electromagnetic, andelastic cases. The most general form of expressing the scattered fieldgenerated by scatterer j due to its exciting field is through theso-called transition matrix operator as follows:

S( r/r _(j) ;a _(j) :{r ₁ , . . . ,r _(j) , . . . ,r _(N) ;a ₁ , . . .,a _(j) , . . . ,a _(N)})= T ( o,m,a) E( r/r _(j) ;a _(j) :{r ₁ , . . .,r _(j) , . . . ,r _(N) ; a ₁ , . . . ,a _(j) , . . . ,a _(N)})   (12)

[0045] In plain words, the transition matrix constitutes a behavioralsignature for the basic individual interaction between the wave and asingle particle suspended in the medium. I call it the particle/mediumsignature. Notice that both the exciting and scattered fields dependupon the spatial/attribute configuration of scatterers, while thetransition operator only depends upon the host medium properties, theparticle j properties, the operating parameters (frequency, temperature,etc) and the nature of the incident physical wave (acoustic,electromagnetic, elastic). It is a representation of the Green'sfunction for a single scatterer in the host medium (single scatteringproblem). When the incident wave is electromagnetic, this inventionemploys the well-known Mie-Scattering model to define the transitionoperator. In the acoustic case the transition operator is obtainedfrom: 1) [5] P. S Epstein and R. R. Carhart, “The Absorption of Sound inSuspensions and Emulsions, I. Water Fog in Air”, J. Acoust. Soc. Am., 25[3] 553-565 (1953), 2) A correction to the thermal diffusion effect dueto surface tension is included from the work of J. C. F. Chow, [14]“Attenuation of Acoustic Waves in Dilute Emulsions and Suspensions”, J.Acoust. Soc. Am., 36 [12] 2395-2401 (1964), 3) A more complete treatmentfor both solid and fluid particles is taken from J. R. Allegra and S. A.Hawley, [19] “Attenuation of Sound in Suspensions and Emulsions: Theoryand Experiments”, J. Acoust. Soc. Am., 51, 1545-1564 (1972), and 4)Viscoelastic behavior of both particles and host medium is modeled alongthe lines described in C. Verdier and M. Piau, [47] “Acoustic WavePropagation in Two-phase viscoelastic fluids: The Case of Polymeremulsions”, J. Acoust. Soc. Am. 101[4], 1868-1876, April 1997 and J. Y.Kim et al, [44] “Dispersive Wave Propagation in a Viscoelastic MatrixReinforced by Elastic Fibers”, J. Acoust. Soc. Am., 95 (3), March 1994.The transition operator for elastic waves is obtained through acombination of the treatment described for the acoustic case and theapproach described in 1) N. G. Einspruch and R. Truell, [7] “Scatteringof a Plane Longitudinal Wave by a Spherical Fluid Obstacle in an ElasticMedium”, J. Acoust. Soc. Am., 32 [2], 214-220, February 1960, and 2) N.G. Einspruch et al, [8] “Scattering of a Plane Transverse Wave by aSherical Obstacle in an Elastic Medium”, J. Appl. Phys., 31 [5],806-817, May 1960.

[0046] In terms of this basic signature and equation 11, the globalbehavior of the ensemble (multiple scattering) will be predicted. Takingthe configurational average on both sides of equation 12 it is obtained:

<S( r/r _(j) ;a _(j))>= T ( o,m,a _(j))<E( r/r _(j) ;a _(j))>  (13)

[0047] Replacing equation 13 in equation 11 the following importantrelation is valid: $\begin{matrix}{{\langle{\underset{\_}{E}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\rangle} = {{\underset{\_}{I}\left( \underset{\_}{r} \right)} + {\int_{Vi}^{\quad}{\int_{D}^{\quad}{{\underset{\_}{\underset{\_}{T}}\left( {\underset{\_}{o},\quad \underset{\_}{m},\quad {\underset{\_}{a}}_{j}} \right)}{\langle{\underset{\_}{E}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{j}};{\underset{\_}{a}}_{j}} \right)}\rangle}{n\left( {{\underset{\_}{r}}_{j},\quad {\underset{\_}{a}}_{j}} \right)}{g\left( {{\underset{\_}{r}}_{j},\quad {\underset{\_}{a}}_{j},\quad {\underset{\_}{r}}_{i},\quad {\underset{\_}{a}}_{i}} \right)}\quad {{\underset{\_}{a}}_{j}}\quad {v_{j}}}}}}} & (14)\end{matrix}$

[0048] Given a fixed spatial/attribute distribution of scatterers, allfields (external and internal to the particles) are supposed to verifythe Helmholtz's vector wave equation:

(∇² +k ²) F ( r:{r ₁ , . . . ,r _(i) , . . . ,r _(N) ;a ₁ , . . . ,a_(i), . . . , a _(N)})=0 ∀Fε{ I,E,S}  (15)

[0049] The symbol ∇² stands for the Laplacian of a vectorial field and kis the wave number for the particular field. The general solution tothis equation may be expressed as a linear combination of threeindependent discrete sets of vector solutions {L _(nm)} (longitudinal),{M _(nm)} (transversal) and {N _(nm)} (transversal) with arbitraryconstants {f_(nm1),f_(nm2),f_(nm3)} as follows [1]: $\begin{matrix}{{\underset{\_}{F}\left( \underset{\_}{r} \right)} = {\sum\limits_{n = 0}^{\infty}\quad {\sum\limits_{m = {- n}}^{m = n}\quad {\left\lbrack {{f_{nm1}{L_{nm}\left( \underset{\_}{r} \right)}} + {f_{nm2}{M_{nm}\left( \underset{\_}{r} \right)}} + {f_{nm3}{N_{nm}\left( \underset{\_}{r} \right)}}} \right\rbrack \quad {\forall{\underset{\_}{F} \in \left\{ {I,\quad E,\quad S} \right\}}}}}}} & (16)\end{matrix}$

[0050] Purely irrotational fields could be expressed exclusively interms of the {L _(nm)} eigenvectors; purely solenoidal fields could beexpressed exclusively in terms of both {M _(nm)} and {N _(nm)}eigenvectors. All three types of basis vector functions are directedrelated to the discrete set of solutions {ψ(r)_(nm)} for thecorresponding scalar Helmholtz wave equation [1]. For particles withshapes randomly deviating from the spherical shape, sphericalcoordinates are appropriate. For fiber-shaped particles, cylindricalcoordinates are the suitable choice. For the sake of continuing thedisclosure, using spherical coordinates (r,θ,φ), the associated discretesets {L _(nm)} {M _(nm)} {N _(nm)} are:

ψ_(nm)( r )=ψ_(nm)(r,θ,φ)=P _(n) ^(m)(cosθ).z _(n)(kr)l ^(mφ) ∀nε{0,1, .. . ,∞};∀mε{−n, . . . 0 . . . n}L _(nm)=∇ψ_(nm) M _(nm) =∇×rψ _(nm) =L_(nm) ×r N _(nm)=(1/k)∇× M _(nm)   (17)

[0051] The {P_(n) ^(m)(cosθ)} set comprises of the well known associatedLegendre Polynomials; and the {z_(n)} set is either the set of sphericalBessel functions {j_(n)} or the set of Hankel functions {h_(n)}. TheBessel set is chosen for those fields that are to be regular at theorigin (incident field and the internal field of a particle); the Hankelfunctional set is used for those fields that are to verify the radiationcondition at infinity (the scattered fields). Superscripts B or H willbe added to {L _(nm)},{M _(nm)},{N _(nm)} according to the case.

[0052] Let the symbol Ndw stand for “Number of different waves”occurring during the interaction between the exciting wave and theparticle. For instance, in the acoustic problem, as described in U.S.Pat. No. 5,121,629, Ndw would be equal to 3 with the value 1corresponding to the compressional waves, value 2 corresponding to thelongitudinal thermal waves, and the value 3 corresponding to the viscoustransversal waves. Adding then a new subscript k=1, Ndw, and definingk_(k) as the corresponding wave number, the following equalities arevalid in spherical coordinates: $\begin{matrix}\begin{matrix}{{L_{nmk}\left( \underset{\_}{r} \right)} = {{\left\{ {\frac{\partial}{\partial r}{z_{n}\left( {k_{k}r} \right)}{P_{n}^{m}\left( {\cos \quad \theta} \right)}^{\quad m\quad \varphi}} \right\} {\underset{\_}{e}}_{r}} +}} \\{{{\left\{ {{z_{n}\left( {k_{k}r} \right)}\frac{\partial}{r{\partial r}}{P_{n}^{m}\left( {\cos \quad \theta} \right)}^{\quad m\quad \varphi}} \right\} {\underset{\_}{e}}_{\theta}} +}} \\{{\left\{ {\frac{\quad m}{r\quad \sin \quad \theta}{z_{n}\left( {k_{k}r} \right)}{P_{n}^{m}\left( {\cos \quad \theta} \right)}^{\quad m\quad \varphi}} \right\} {\underset{\_}{e}}_{\varphi}}} \\{{M_{nmk}\left( \underset{\_}{r} \right)} = {{\left\{ {\frac{\quad m}{\sin \quad \theta}{z_{n}\left( {k_{k}r} \right)}{P_{n}^{m}\left( {\cos \quad \theta} \right)}^{\quad m\quad \varphi}} \right\} {\underset{\_}{e}}_{\theta}} +}} \\{{\left\{ {{z_{n}\left( {k_{k}r} \right)}\frac{\partial}{\partial\theta}{P_{n}^{m}\left( {\cos \quad \theta} \right)}^{\quad m\quad \varphi}} \right\} {\underset{\_}{e}}_{\varphi}}} \\{{N_{nmk}\left( \underset{\_}{r} \right)} = {{\left\{ {\frac{n\left( {n + 1} \right)}{k_{k}r}{z_{n}\left( {k_{k}r} \right)}{P_{n}^{m}\left( {\cos \quad \theta} \right)}^{\quad m\quad \varphi}} \right\} {\underset{\_}{e}}_{r}} +}} \\{{{\left\{ {\frac{\partial}{k_{k}r{\partial r}}{{rz}_{n}\left( {k_{k}r} \right)}\frac{\partial}{\partial\theta}{P_{n}^{m}\left( {\cos \quad \theta} \right)}^{\quad m\quad \varphi}} \right\} {\underset{\_}{e}}_{\theta}} +}} \\{{\left\{ {\frac{\partial}{k_{k}r{\partial r}}{{rz}_{n}\left( {k_{k}r} \right)}\frac{\quad m}{\sin \quad \theta}{P_{n}^{m}\left( {\cos \quad \theta} \right)}^{\quad m\quad \varphi}} \right\} {\underset{\_}{e}}_{\varphi}}}\end{matrix} & (18)\end{matrix}$

[0053] It is worth noting that equations 5 and 6 of prior art U.S. Pat.No. 5,121,629 are special versions of our generic equation 16 inspherical coordinates (upon use of equation 18) for the simpler acousticcase where the incident wave is purely longitudinal (a field ofgradients). Let us now apply our equation 16 to the exciting fieldsurrounding particle i located at position r _(i) as well as to thefield scattered by particle j located at position r _(j), and takeconfigurational averages: $\begin{matrix}\begin{matrix}{{\langle{\underset{\_}{E}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{i}};{\underset{\_}{a}}_{i}} \right)}\rangle} = {\sum\limits_{n = 0}^{\infty}{\sum\limits_{m = {- n}}^{m = n}{\sum\limits_{k = 1}^{Ndw}\left\lbrack {{{\langle{e_{{nmk}\quad 1}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{L_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} +} \right.}}}} \\{{{{\langle{e_{{nmk}\quad 2}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{M_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} +}} \\\left. {{\langle{e_{{nmk}\quad 3}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{N_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} \right\rbrack \\{{\forall{\underset{\_}{r}:{{{\underset{\_}{r} - {\underset{\_}{r}}_{i}}} \leq b}}}}\end{matrix} & (19) \\\begin{matrix}{{\langle{\underset{\_}{S}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{j}};{\underset{\_}{a}}_{j}} \right)}\rangle} = {\sum\limits_{n = 0}^{\infty}{\sum\limits_{m = {- n}}^{m = n}{\sum\limits_{k = 1}^{Ndw}\left\lbrack {{{\langle{s_{{nmk}\quad 1}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}{L_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)}} +} \right.}}}} \\{{{{\langle{s_{{nmk}\quad 2}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}{M_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)}} +}} \\\left. {{\langle{s_{{nmk}\quad 3}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}{N_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)}} \right\rbrack \\{{\forall{\underset{\_}{r}:{{{\underset{\_}{r} - {\underset{\_}{r}}_{j}}} \geq b}}}}\end{matrix} & (20)\end{matrix}$

[0054] The symbol b stands for the distance of closest approach betweenparticles i and j (b≧(d_(j)+d_(i))/2). Following the expansions inequations 19 and 20, applying a suitable set of boundary conditionsaround particle j, and adding the polarization index q=1,2,3, theabstract nature of the operator T(o,m,a _(j)) materializes into amulti-dimensional infinite matrix that relates each scatteringcoefficient in equation 20 with all the coefficients in equation 19 as:$\begin{matrix}\begin{matrix}{{\langle{s_{nmkq}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle} = {\sum\limits_{r = 0}^{\infty}{\sum\limits_{s = {- r}}^{s = r}{\sum\limits_{t = 1}^{Ndw}{\sum\limits_{u = 1}^{3}{{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{rstu}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}}}}}} \\{= {\sum\limits_{rstu}{{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{rstu}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}}}\end{matrix} & (21)\end{matrix}$

[0055] In the electromagnetic case, due to the purely solenoidalcharacter of the wave and the spherical symmetry involved, only T_(n112)^(n112) (o,m,p) and T_(n113) ^(n113) (o,m,p _(j)) are needed and theycoincide with the well known Mie coefficients an and b, in the prior art[27]. In the acoustic and elastic cases both longitudinal andtransversal waves exist. Applying the single summation convention usedin equation 21, equation 20 now becomes: $\begin{matrix}\begin{matrix}{{\langle{\underset{\_}{S}\left( {{\underset{\_}{r}/{\underset{\_}{r}}_{j}};{\underset{\_}{a}}_{j}} \right)}\rangle} = {\sum\limits_{nmk}\left\lbrack \left\{ {\sum\limits_{rst}{{T_{{rst}\quad 1}^{{nmk}\quad 1}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{{rst}\quad 1}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}} \right\} \right.}} \\{{{L_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)} + {\sum\limits_{rst}{T_{{rst}\quad 2}^{{nmk}\quad 2}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}}}} \\{{\left. {\langle{e_{{rst}\quad 2}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle} \right\} {M_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)}} +} \\{\left\{ {\sum\limits_{rst}{T_{{rst}\quad 3}^{{nmk}\quad 3}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right){\langle{e_{{rst}\quad 3}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}} \right\}} \\\left. {N_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)} \right\rbrack\end{matrix} & (22)\end{matrix}$

[0056] And, consequently, equation 14 converts to: $\begin{matrix}{{\sum\limits_{nmk}\left\lbrack {{{\langle{e_{{nmk}\quad 1}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{L_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {{\langle{e_{{nmk}\quad 2}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{M_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {{\langle{e_{{nmk}\quad 3}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{N_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}} \right\rbrack} = {{\sum\limits_{nmk}\left\lbrack {{{i_{{nmk}\quad 1}\left( {\underset{\_}{r}}_{i} \right)}{L_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {{i_{{nmk}\quad 2}\left( {\underset{\_}{r}}_{i} \right)}{M_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {{i_{{nmk}\quad 3}\left( {\underset{\_}{r}}_{i} \right)}{N_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}} \right\rbrack} + {\int_{Vi}{\int_{Di}{{n\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j}} \right)}{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i},{\underset{\_}{a}}_{i}} \right)}{\sum\limits_{nmk}{\left\lbrack {{\left\{ {\sum\limits_{{rst}\quad 1}{{T_{{rst}\quad 1}^{{nmk}\quad 1}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{{rst}\quad 1}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}} \right\} {L_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)}} + {\left\{ {\sum\limits_{rst}{{T_{{rst}\quad 2}^{{nmk}\quad 2}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{{rst}\quad 2}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}} \right\} {M_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)}} + {\left\{ {\sum\limits_{rst}{{T_{{rst}\quad 3}^{{nmk}\quad 3}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{{rst}\quad 3}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}} \right\} {N_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)}}} \right\rbrack {{\underset{\_}{a}}_{j}}{v_{j}}}}}}}}} & (23)\end{matrix}$

[0057] In the above equation the eigenvectors are expanded around thelocation of particle i in both the exciting field and the incidentfield. However, in the terms associated with the scattered fields, theeigenvectors are expanded around the location of particle j. It is worthnoting as well that the latter eigenvectors are expressed in terms ofHankel functions while the former are expressed in terms of Besselfunctions. In order to take subsequent advantage of the orthogonality ofthese vector functions, all of them need to be expressed in terms of thesame type of functions and around a single location. For the purpose, Iuse the translational addition theorem for spherical vector wavefunctions well known in the prior literature which can be stated as:$\begin{matrix}\begin{matrix}{{L_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)} = {\sum\limits_{v = 0}^{\infty}{\sum\limits_{\mu = {- v}}^{v}{{A_{v\quad \mu}^{{nmk}\quad 1}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{L_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)}}}}} \\{= {\sum\limits_{v\quad \mu}{{A_{v\quad \mu}^{{nmk}\quad 1}\left( {{\underset{\_}{r}}_{j} - \underset{\_}{r_{i}}} \right)}{L_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}}} \\{{M_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)} = {{\sum\limits_{v\quad \mu}{{A_{v\quad \mu}^{{nmk}\quad 2}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{M_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}} +}} \\\left. {{A_{v\quad \mu}^{{nmk}\quad 3}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{N_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} \right\rbrack \\{{N_{nmk}^{H}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{j}} \right)} = {{\sum\limits_{v\quad \mu}{{A_{v\quad \mu}^{{nmk}\quad 3}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{M_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}} +}} \\\left. {{A_{v\quad \mu}^{{nmk}\quad 2}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{N_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} \right\rbrack\end{matrix} & (24)\end{matrix}$

[0058] Where the translational coefficients A_(vμ) ^(nmkq) are given by:$\begin{matrix}{A_{v\quad \mu}^{nmkq} = {\sum\limits_{p}{{a_{q}\left( {n,m,v,\mu,p} \right)}{h_{p}\left( {k_{k}r_{ij}} \right)}{P_{p}^{m - \mu}\left( {\cos \quad \theta_{ij}} \right)}{\exp\left( {{\left( {m - \mu} \right)}\varphi_{ij}} \right.}}}} & (25)\end{matrix}$

[0059] The index p varies over the range |n−v|,|n−v+2, . . . n+v and thecoefficients a_(q)(n,m,v,μ,p) are:a₁(n, m, v, μ, p) = (−1)^(μ)^(v + p − n)(2v + 1)a(m, n−μ, vp)a₂(n, m, v, μ, p) = (−1)^(μ)^(v + p − n)a(n, v, p)a(m, n−μ, vp)a₃(n, m, v, μ, p) = (−1)^(μ)^(v + p − n)b(n, v, p)a(m, n−μ, vp, p − 1)a(n, v, p), b(n, v, p), a(m, n, −μ, vp)  anda(m, n, −μ, v  p, p − 1)  are:${a\left( {n,v,p} \right)} = {{{\left\lbrack {{2{v\left( {v + 1} \right)}\left( {{2v} + 1} \right)} + {\left( {v + 1} \right)\left( {n - v + p + 1} \right)\left( {n + v - p} \right)} - {{v\left( {v - n + p + 1} \right)}\left( {n + v + p + 2} \right)}} \right\rbrack/\left\lbrack {2{v\left( {v + 1} \right)}} \right\rbrack}{b\left( {n,v,p} \right)}} = {{\left\lbrack {\left( {n + v + p + 1} \right)\left( {v - n + p} \right)\left( {n - v + p} \right)\left( {n + v - p + 1} \right)} \right\rbrack^{1/2}{\left( {{2v} + 1} \right)/\left\lbrack {2{v\left( {v + 1} \right)}} \right\rbrack}{a\left( {m,{n{{,{- \mu},v}}p}} \right)}} = {{\left( {- 1} \right)^{m + \mu}{\left( {{2p} + 1} \right)\left\lbrack \frac{{\left( {n + m} \right)!}{\left( {v + \mu} \right)!}{\left( {p - m - \mu} \right)!}}{{\left( {n - m} \right)!}{\left( {v - \mu} \right)!}{\left( {p + m + \mu} \right)!}} \right\rbrack}^{1/2}\left( \frac{{n\ldots}\quad {v\ldots}\quad p}{{{m\ldots}\quad {\mu \ldots}} - \left( {m + \mu} \right)} \right)\left( \frac{{n\ldots}\quad {v\ldots}\quad p}{0\ldots \quad 0\ldots \quad 0} \right){a\left( {m,{n{{,{- \mu},v}}p},q} \right)}} = {\left( {- 1} \right)^{m + \mu}{\left( {{2p} + 1} \right)\left\lbrack \frac{{\left( {n + m} \right)!}{\left( {v + \mu} \right)!}{\left( {p - m - \mu} \right)!}}{{\left( {n - m} \right)!}{\left( {v - \mu} \right)!}{\left( {p + m + \mu} \right)!}} \right\rbrack}^{1/2}\left( \frac{{n\ldots}\quad {v\ldots}\quad p}{{{m\ldots}\quad {\mu \ldots}} - \left( {m + \mu} \right)} \right)\left( \frac{{n\ldots}\quad {v\ldots}\quad p}{0\ldots \quad 0\ldots \quad 0} \right)}}}}$

[0060] The symbol$\left( \frac{j_{1}\ldots \quad j_{2}\ldots \quad j_{3}}{{m_{1}\ldots \quad m_{2}\ldots} - \left( {m_{1} + m_{2}} \right)} \right)$

[0061] is the well-known Wigner 3j symbol [6]. Upon application ofequations 24, equation 23 transforms to: $\begin{matrix}{{\sum\limits_{nmk}\left\lbrack {{{\langle{e_{{nmk}\quad 1}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{L_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {{\langle{e_{{nmk}\quad 2}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{M_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {{\langle{e_{{nmk}\quad 3}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle}{N_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}} \right\rbrack} = {{\sum\limits_{nmk}\left\lbrack {{{i_{{nmk}\quad 1}\left( {\underset{\_}{r}}_{i} \right)}{L_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {{i_{{nmk}\quad 2}\left( {\underset{\_}{r}}_{i} \right)}{M_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {{i_{{nmk}\quad 3}\left( {\underset{\_}{r}}_{i} \right)}{N_{nmk}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}} \right\rbrack} + {\int_{Vi}{\int_{D}{{n\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j}} \right)}{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i},{\underset{\_}{a}}_{i}} \right)}{\sum\limits_{nmk}{\left\lbrack {{\left\{ {\sum\limits_{rst}{{T_{{rst}\quad 1}^{{nmk}\quad 1}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{{rst}\quad 1}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}} \right\} {\sum\limits_{v\quad \mu}{A_{v\quad \mu}^{{nmk}\quad 1}{L_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}}} + {\left\{ {\sum\limits_{rst}{{T_{{rst}\quad 2}^{{nmk}\quad 2}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{{rst}\quad 2}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}} \right\} {\sum\limits_{v\quad \mu}\left\lbrack {{A_{v\quad \mu}^{{nmk}\quad 2}{M_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {A_{v\quad \mu}^{{nmk}\quad 3}{N_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}} \right\rbrack}} + {\left\{ {\sum\limits_{rst}{{T_{{rst}\quad 3}^{{nmk}\quad 3}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{p}}_{j}} \right)}{\langle{e_{{rst}\quad 3}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}}} \right\} {\sum\limits_{v\quad \mu}\left\lbrack {{A_{v\quad \mu}^{{nmk}\quad 3}{M_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}} + {A_{v\quad \mu}^{{nmk}\quad 2}{N_{v\quad \mu \quad k}^{B}\left( {\underset{\_}{r} - {\underset{\_}{r}}_{i}} \right)}}} \right\rbrack}}} \right\rbrack {{\underset{\_}{a}}_{j}}{v_{j}}}}}}}}} & (26)\end{matrix}$

[0062] It is now straightforward, applying the mutual orthogonality ofthe functional sets {L _(nmk)},{M _(nmk)},{N _(nmk)}, to arrive to thefollowing equation: $\begin{matrix}{{{\langle{e_{nmkq}\left( {{\underset{\_}{r}}_{i};{\underset{\_}{a}}_{i}} \right)}\rangle} = {{i_{nmkq}\left( {\underset{\_}{r}}_{i} \right)} + {\sum\limits_{rstu}{\int_{Vi}{{n\left( {\underset{\_}{r}}_{j} \right)}\left\{ {\int_{D}{{f_{n}\left( {\underset{\_}{a}}_{j} \right)}{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i},{\underset{\_}{a}}_{i}} \right)}{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{\langle{e_{rstu}\left( {{\underset{\_}{r}}_{j};{\underset{\_}{a}}_{j}} \right)}\rangle}{{\underset{\_}{a}}_{j}}}} \right\} {\sum\limits_{v\quad \mu}{{A_{v\quad \mu}^{nmkq}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{v_{j}}}}}}}}}} & (27)\end{matrix}$

[0063] Taking now the configurational average with respect to theattributes of particle i I find: $\begin{matrix}{{\langle{\langle{e_{nmkq}\left( {\underset{\_}{r}}_{i} \right)}\rangle}\rangle} = {{{i_{nmkq}\left( {\underset{\_}{r}}_{i} \right)} + {\sum\limits_{rstu}{\int_{Vi}{{n\left( {\underset{\_}{r}}_{j} \right)}{\langle{T_{rdtu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{r}}_{j}} \right)}\rangle}{\langle{\langle{e_{rstu}\left( {\underset{\_}{r}}_{j} \right)}\rangle}\rangle}{\sum\limits_{v\quad \mu}{{A_{v\quad \mu}^{nmkq}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{v_{j}}{\langle{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{r}}_{j}} \right)}\rangle}}}}}}} = {{\int_{D2}{{f_{n}\left( {\underset{\_}{a}}_{j} \right)}\left( {g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i}} \right)}\rangle \right.{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{{\underset{\_}{a}}_{j}}{\langle{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i}} \right)}\rangle}}} = {\int_{D1}{{f_{n}\left( {\underset{\_}{a}}_{i} \right)}{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i},{\underset{\_}{a}}_{i}} \right)}{{\underset{\_}{a}}_{i}}}}}}} & (28)\end{matrix}$

[0064] I have applied the equalities <<e_(rstu)(r _(j);a_(j))>>=<<e_(rstu)(r _(j);a _(i))>>=<<e_(rstu)(r _(j))>> that resultfrom the non-segregation hypothesis. Domains D1 and D2 are governed bythe non-interpenetrability of the particles because some properties(like size and shape) relate to the space occupancy of the particles. Itis worth noting the radical difference between equation 28 (for the casewhen only the attribute size is considered) and the prior art, e.g.,[36] Varadan et al in page 379 of “A Multiple Scattering Theory forElastic Wave Propagation in Discrete Random Media”, J. Acoust. Soc. Am.,77(2), February 1985, where the authors suggest that a distribution ofsizes could be dealt with by replacing the transition operator T by⟨T⟩ = ∫₀^(∞)T(d)f_(n)(d)d.

[0065] Even recalculating the correlation integral for each size of thedistribution as the authors suggest as more accurate, it is clear thatsuch an approach can only constitute a coarse approximation to thecomplex intricate relation between <T>,T,f_(n),<g>, and g shown inequation 28. They agree only when <g(r _(j),a _(j),r _(i),a _(i))>≡1--acondition not true even under the “hole correction”, i.e. when CF(r)≡1.The main reason for this disparity resides in the prior-art neverrecognizing that when the population has a distribution of attributes,so it does the arbitrary particle chosen as fixed in space so as todevelop the stochastic field equations. Consequently, the secondconfigurational average with respect to the attributes is paramount toarrive to the correct final dispersion equations.

[0066] Before proceeding towards the final dispersion equations for theparticulate, let us see how equation 28 transforms for simpler casesalready contemplated in the prior art. For the simpler but importantbasic case in which all scatterers are identical with attribute vector a₀, i.e.—using the Dirac's δ function—f_(n)(a _(k))=δ(a _(k)−a ₀), thefollowing equalities are valid: $\begin{matrix}\begin{matrix}{{\langle{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i}} \right)}\rangle} = {\int_{D1}{{\delta \left( {{\underset{\_}{a}}_{i} - {\underset{\_}{a}}_{0}} \right)}{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i},{\underset{\_}{a}}_{i}} \right)}{{\underset{\_}{a}}_{i}}}}} \\{= {g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i},{\underset{\_}{a}}_{0}} \right)}} \\{{\langle{T_{rstq}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{r}}_{j}} \right)}\rangle} = {\int_{D2}{{\delta \left( {{\underset{\_}{a}}_{j} - {\underset{\_}{a}}_{0}} \right)}{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{j},{\underset{\_}{r}}_{i},{\underset{\_}{a}}_{0}} \right)}{T_{rstq}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{j}} \right)}{{\underset{\_}{a}}_{j}}}}} \\{= {{g\left( {{\underset{\_}{r}}_{j},{\underset{\_}{a}}_{0},{\underset{\_}{r}}_{i},{\underset{\_}{a}}_{0}} \right)}{T_{rstq}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},a_{0}} \right)}}}\end{matrix} & (29)\end{matrix}$

[0067] Further assuming the particle concentration is spatially uniform,i.e. n(r _(j))=n₀, and the pair-correlation being a function of only thedistance between the identical particles and their identical size D,equation 28 becomes: $\begin{matrix}{{\langle{{\langle{e_{nmkq}\left( {\underset{\_}{r}}_{i} \right)}\rangle}}\rangle} = {{i_{nmkq}\left( {\underset{\_}{r}}_{i} \right)} + {n_{0}{\sum\limits_{rstq}{{T_{rstq}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{a}}_{0}} \right)}{\int_{V_{1}}{{g\left( {{{{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}}},D} \right)}{\langle{{\langle{e_{rstq}\left( {\underset{\_}{r}}_{j} \right)}\rangle}}\rangle}{\sum\limits_{v\quad \mu}{{A_{v\quad \mu}^{nmkq}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{v_{j}}}}}}}}}}} & (30)\end{matrix}$

[0068] The above equation is found in the prior-art literature foridentical scatterers with the following caveats: a) the configurationalaveraging with respect to a _(i) (second angular bracket) does notappear though, being a _(i)=a _(j)∀i,j when scatterers are identical, itis immaterial; b) the concept behind the existence of subscript k (to,e.g., contemplate thermal and viscous-inertial waves) is definitelyinexistent in the prior art. This latter feature is essential for anaccurate modeling of multiple scattering in the acoustical andelastodynamic cases. It is clear then that the description of themultiple scattering phenomena in particulate composites provided by thisinvention is considerably more general and accurate than those in theexisting prior art. Besides, as it should be, upon the consideration ofsimpler specific cases, our equations reduce to those well known by theskilled in the art. The assertion regarding accuracy will be furtherproved with experimental data.

[0069] Equation 28 constitutes a linear system of equations with<<e_(nmkq)(r _(i))>> as the unknowns. In order to solve it, I need tomake some practical assumptions. I start assuming that the coherentexciting field <<e_(nmkq)(r _(i))>> is that of a plane wave with aneffective propagation constant K_(k). In symbols: I look for a solutionof the form <<e_(nmkq)(r _(i))>>=X_(nmkq)exp(iK_(k)z_(i)). Replacing inequation 28 it is obtained: $\begin{matrix}{{X_{nmkq}{\exp \left( {\quad K_{k}z_{i}} \right)}} = {{i_{nmkq}\left( {\underset{\_}{r}}_{i} \right)} + {n_{0}{\sum\limits_{rstu}{X_{rstu}{\exp \left( {\quad K_{t}z_{i}} \right)}{\int_{Vi}{{\langle{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{r}}_{j}} \right)}\rangle}{\exp \left( {\quad {K_{t}\left( {z_{j} - z_{i}} \right)}} \right)}{\sum\limits_{v\quad \mu}{{A_{v\quad \mu}^{nmkq}\left( {{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}} \right)}{v_{j}}}}}}}}}}} & (31)\end{matrix}$

[0070] Using equation 25 for the Coefficients A_(vμ) ^(nmkq) (r _(j)−r_(i)) and noting the axial symmetry of the problem that renders onlyterms independent of the azimuthal angle (m=μ) as contributors to theintegral, it is obtained: $\begin{matrix}{{{X_{nmkq}{\exp \left( {\quad K_{k}z_{t}} \right)}} = {{i_{nmkq}\left( {\underset{\_}{r}}_{i} \right)} + {n_{0}{\sum\limits_{rstu}{X_{rstu}{\exp \left( {\quad K_{t}z_{t}} \right)}{\sum\limits_{v\quad \mu}{\sum\limits_{p}{{a_{q}\left( {n,m,v,\mu,p} \right)}{I_{s\quad \mu \quad p}\left( {K_{t},k_{t}} \right)}}}}}}}}}{{I_{s\quad \mu \quad p}\left( {K_{t},k_{t}} \right)} = {\delta_{s\quad \mu}{\int_{Vi}{{\langle{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{r}}_{j}} \right)}\rangle}{\exp \left( {\quad {K_{t}\left( {z_{j} - z_{i}} \right)}} \right)}{h_{p}\left( {k_{t}{{{\underset{\_}{r}}_{j} - {\underset{\_}{r}}_{i}}}} \right)}{P_{p}\left( {\cos \quad \theta_{ij}} \right)}{\exp \left( {{\left( {s - \mu} \right)}\varphi_{ij}} \right)}{v_{j}}}}}}} & (32)\end{matrix}$

[0071] Where δ_(ij) is the Kronecker symbol which value is 1 for i=j and0 for i≠j. Because both the unknown coherent wave in the composite andthe set of basis functions verify the Helmholtz equation, the integralcan be rewritten as: $\begin{matrix}{{I_{s\quad \mu \quad p}\left( {K_{t},k_{t}} \right)} = {\delta_{s\quad \mu}{\int_{r \geq b}{{\langle{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{r}}_{j}} \right)}\rangle}{\exp \left( {\quad K_{t}z} \right)}{h_{p}\left( {k_{t}r} \right)}{P_{p}\left( {\cos \quad \theta_{ij}} \right)}{v}}}}} \\{= {\left\{ {1/\left( {K_{t}^{2} - k_{t}} \right)} \right\} {\int_{r \geq b}{{{\langle{T_{rstu}^{nmkq}\left( {\underset{\_}{o},\underset{\_}{m},{\underset{\_}{r}}_{j}} \right)}\rangle}\left\lbrack {{{\exp \left( {\quad K_{t}z} \right)}{\nabla^{2}\left\{ {{h_{p}\left( {k_{t}r} \right)}{P_{p}\left( {\cos \quad \theta_{ij}} \right)}} \right\}}} - {{\nabla^{2}\left\{ {\exp \left( {\quad K_{t}z} \right)} \right\}}{h_{p}\left( {k_{t}r} \right)}{P_{p}\left( {\cos \quad \theta_{ij}} \right)}}} \right\rbrack}{v}}}}}\end{matrix}$

[0072] And using Gauss's theorem of divergence so as to convert thevolume integral into a surface integral, the value for the integral canbe expressed as the sum of three terms I_(p) ^(c), I_(p) ^(d) and I_(p)^(e) as follows: $\begin{matrix}\begin{matrix}{{I_{p}^{c}\left( {K_{t},k_{t}} \right)} = {\left\{ \frac{2\pi \quad i}{k_{t}^{2}\left( {K_{t} - k_{t}} \right)} \right\} i^{p}{\exp \left( {\quad k_{t}z_{i}} \right)}{\exp \left( {{- }\quad K_{t}z_{i}} \right)}}} \\{{I_{p}^{d}\left( {K_{t},k_{t},b} \right)} = {{- \left\{ \frac{4\pi \quad b^{2}}{\left( {K_{t}^{2} - k_{t}^{2}} \right)} \right\}}{i^{p}\left\lbrack {{\frac{\partial}{\partial b}{h_{p}\left( {k_{t}b} \right)}{j_{p}\left( {K_{t}b} \right)}} -} \right.}}} \\\left. {{h_{p}({kb})}\frac{\partial}{\partial b}{j_{p}\left( {K_{t}b} \right)}} \right\rbrack \\{{I_{p}^{e}\left( {K_{t},k_{t},b} \right)} = {\int_{b}^{\infty}{{r^{2}\left\lbrack {{\langle{T_{rstq}\left( {\underset{\_}{o},\underset{\_}{m},r} \right)}\rangle} - 1} \right\rbrack}{h_{p}\left( {k_{t}r} \right)}{j_{p}\left( {K_{t}r} \right)}{r}}}}\end{matrix} & (33)\end{matrix}$

[0073] It is worth noting once more the difference in this case betweenour I_(p) ^(e)(K_(t),k_(t),b)and the corresponding term in theprior-art, e.g. term M_(p) (k,K|b) in page 554 of [28] L. Tsang et al,“Multiple Scattering of Acoustic Waves by Random Distribution ofDiscrete Spherical Scatterers with the Quasicrystalline andPercus-Yevick Approximation”, J. Acoust. Soc. Am. 71(3), March 1982. Inthis prior-art publication our term <T_(rstq)(o,m,r)> does not exit, itis simply replaced by the transition matrix times the pair correlationfunction, a replacement that is only valid for identical scatterers.

[0074] Finally, without loss of generality, the coordinate system can beselected such that the z-axis is in the direction of propagation of theincident wave and the y-axis is in the direction of the displacementvector for the transversal component of the incident wave. Consequently,a generic representation for a plane incident wave (acoustic,electromagnetic or elastic) of unity intensity with respect to the pointr _(i) is:

I ( r−r _(i))=exp(ik ₂ z _(i))exp(ik ₂(z−z _(i))) e _(y)+exp(ik ₁ z_(i))exp(ik ₁(z−z _(i))) e_(z) r=xe _(x) +ye _(y) +ze_(z)   (34)

[0075] Let us now explicitly express i_(nmkq)(r _(i)) in sphericalcoordinates [1].

i _(nmlq)( r _(i))={exp(ik ₁ z _(i))/ik ₁}(2n+1)i ^(n)δ_(m0)δ_(q1)

i _(nm21)( r _(i)=0 i _(nm22)( r _(i))=({exp(ik ₂ z _(i))/2}{(2n+1)i^(n)/(n(n+1)}{δ_(m1) −n(n+1)δ_(m,−1)}

i _(nm23)( r _(i))={exp(ik ₂ z _(i))/2}{(2n+1)i ^(n)/(n(n+1)}{δ_(m1)+n(n+1)δ_(m,−1)}  (35)

[0076] When placed into equation 32 the first term in equation 33 leadsto the cancellation of the incident wave as per the well-knownEwald-Oseen extinction theorem expressed in the following equation:$\begin{matrix}{0 = {{i_{nmkq}\left( {\underset{\_}{r}}_{i} \right)} + {n_{0}{\sum\limits_{rstu}{X_{rstu}{\sum\limits_{v\quad \mu}{\sum\limits_{p}{{a_{q}\left( {n,m,v,\mu,p} \right)}\left\{ \frac{2\pi \quad i}{k_{t}^{2}\left( {K_{t} - k_{t}} \right)} \right\} i^{p}}}}}}}}} & (36)\end{matrix}$

[0077] From the above equation the multiple-scattering propagationconstants for the particulate composite are obtained as a function ofthe propagation constants of the host medium, the particleconcentration, and the amplitudes X_(rstu). Placing now the second andthird terms of equation 33 into equation 32, the following genericexpression for the Lorentz-Lorenz law is obtained: $\begin{matrix}{X_{nmkq} = {{i_{nmkq}\left( {\underset{\_}{r}}_{t} \right)} + {n_{0}{\sum\limits_{rstu}{X_{rstu}{\sum\limits_{v\quad \mu}{\sum\limits_{p}{{a_{q}\left( {n,m,v,\mu,p} \right)}\left\{ {{I_{p}^{d}\left( {K_{t},k_{t}} \right)} + {I_{p}^{e}\left( {K_{t},k_{t}} \right)}} \right\}}}}}}}}} & (37)\end{matrix}$

[0078] The last equation constitutes an infinite system of simultaneoushomogeneous equations in complex numbers for the unknown amplitudesX_(nmkq) that can be truncated to a finite system upon employment of asuitable convergence criterion. A non-trivial solution will exist onlywhen the determinant of the system is zero. The propagation constantsfor the particulate system are then determined by finding the value ofK, for which the determinant is null. Any suitable algorithm foraccurately finding the complex root of a complex function could be usedin this invention (e.g. from “Collected Algorithms from ACM”, TheAssociation for Computing Machinery, or W. H. Press et al, “NumericalRecipes, The Art of Scientific Computing”, Cambridge Press, 1988).Knowing the propagation constants for the composite, all its attributes(attenuation, scattering, velocity, physical properties, etc) can bestraightforwardly determined. For instance, the expected value (χ) forthe adiabatic compressibility of the composite can be calculated fromthe formula${\langle\chi\rangle} = \frac{K^{2}\chi_{m}\rho_{m}}{k_{m}^{2}{\langle\rho\rangle}}$

[0079] where subscript m stands for host medium, <ρ> is the averagedensity of the composite, and K and k are the composite and host mediumpropagation constants respectively.

[0080] Experimental Validation of the METAMODEL™ Prediction Engine 21

[0081] Well-controlled experimental data in concentrated particulatesare not abundant in the technical literature. Some data forelectromagnetic and elastic waves, and their successful comparison withpredictions from specific versions of the analytical theory can be foundin [31,32,33,36,37,42]. Because the present prediction engine is theresult of a generic employment of the analytical theory plus the novelfeatures already described, those prior successes can be readilyconsidered as a partial validation for the prediction engine of thisinvention. However, this by itself would not be —obviously—enough. Themost appropriate field to fully validate the current invention with itsnovel features is in Acoustics because of: a) the diversity of physicalphenomena occurring when a sound wave interacts with particulates, b)accurate acoustic data have been made available only during the lastdecade due to the recent commercial advent of acoustic spectrometers,and c) it is the successful integration of viscous-inertial, thermaldiffusion, surface tension, and viscoelasticity phenomena into thefundamental statistical treatment of multiple-scattering whatconstitutes the gist for the prediction engine of the current invention.Consequently, I disclose herein a considerable amount of acousticexperimental data that solidly validate the operability of theMETAMODEL™ Prediction Engine.

[0082] Acoustical data in suspensions and emulsions were gathered usingthe acoustic spectrometer function of a particle size analyzer(ULTRASIZER™) manufactured by Malvern Instruments Ltd of the UnitedKingdom. In this instrument, the attenuation coefficient of ultrasoundwaves in a frequency range from 1 MHz to 150 MHz is measured withaccuracy better than 5%. Concentration ladders from low concentrations(C_(v)≦1%), where both the standard single-scattering model and theprediction engine of the current invention perfectly agree, up toconcentrations as high as 30% by volume were carefully made avoidingaggregation. For every measured sample, only measurements forconcentrations above 10% are shown in this disclosure because—in theacoustics field—only at such high concentrations is when the drasticimprovement on accuracy achieved by the current invention is easilyvisible.

[0083]FIG. 3 and FIG. 4 correspond to an extremely fine dispersion ofsilica particles (KLEBOSOL® 30R12) used in the semiconductor industryfor silicon wafer polishing. The manufacturer gives only an approximatevalue for the mean size of around 12 nm. Using the ULTRASIZER™ at theoptimal concentration for maximum accuracy, its actual mean size wasdetermined as 19 nm. With this value for the size, the attenuationspectrum at high concentrations was calculated with our predictionengine and compared with the experimental spectra. Predictions based onthe standard single-scattering model are also shown on the graphs inorder to assess the significance of multiple-scattering effects. FIG. 3and FIG. 4 correspond to C_(v)=11.4% and C_(v)=17.1% respectively.Agreement between experimental data and METAMODEL™ predictions is within10% throughout the frequency range, with a typical accuracy of about 5%.Nissan Chemical America Corp. manufactures silica dispersions with verynarrow size distribution and well-controlled mean size. FIG. 5 and FIG.6 correspond to a product with a size specified as 100 nm±30 nm. Theactual size was determined with the ULTRASIZER™ at a concentration formaximum accuracy obtaining 105 nm. With this value for the size, theprediction engine for this invention was run at high concentrations andcompared with the actual measured spectra. FIG. 5 and FIG. 6 are for 10%and 30% respectively. Agreement is again excellent (well below 10%).FIG. 7 and FIG. 8 are for a product specified by the manufacturer as 300nm±30 nm. The ULTRASIZER™, operated at the optimum concentration,delivered 310 nm. Using this value for the size, FIG. 7 and FIG. 8 showa very good agreement between METAMODEL™ and the actual spectra taken at16% and 24%. The last Nissan product is stated to have a mean size of450 nm±30 nm . FIG. 9 and FIG. 10 are, by now, self-explanatory,confirming once more the predicting power of the current invention.Finally, FIG. 11 corresponds to a broad distribution of soy oil dropletsin water at about 11% concentration. As expected for emulsions, multiplescattering effects at this relatively high concentration are still notvery significant [48]. Notwithstanding, there is a visible smallsystematic deviation of the actual measured spectrum from thesingle-scattering prediction at the low frequency end. The METAMODEL™Prediction Engine once more demonstrates its prediction power byfollowing almost perfectly (better than 5%) the experimental data at lowfrequency.

[0084] Detailed Description/Operation of the SCATTERER™ GraphicalInterface 23

[0085] Due to the remarkable advance of computer technology during thelast two decades, de facto industry standards for graphical userinterfaces have developed and are well known and employed in theimplementation of any computer-based system. Hierarchies of windows,icons, option and command buttons, checkbox controls, context menus,image controls, tab-strip controls, tables, graphs, tool-tips, on-linehelp, etc are common place in the attempt to improve the ergonomics ofhuman interaction with machines. Upon installation of the attached CD-Rin a suitable computer, the “look and feel” of a preferred embodimentfor this invention can be experienced. Virtually an infinite number ofinterfaces with different “looks and feels” could be designed, and stillall of them achieving the essential object of this important componentfor the current invention. With that proviso, and for the sake of abetter disclosure of the features and operation of this invention, Iwill describe the structure, materials, and operation of the GraphicalInterface component 23 by direct reference to the preferred embodiment.

[0086] Interaction with the human operator starts by allowing her toselect amid three suites: Acoustics, Electromagnetics, andElastodynamics. After the desired selection, the appearance of thesesuites to the user is basically the same, except for some differentvariables and default conditions corresponding to well-establishedpractices in each field. For instance, in acoustics it is customary torefer to the frequency of the wave, while in electromagnetics it is thewavelength the used equivalent operating attribute; in acoustics thefrequency is most commonly varied and the scattering angle is fixed,while in light scattering exactly the opposite is normally practiced,and so forth. The concept of a tool-tip-text (brief explanation of thepurpose of each control) is used throughout the interface to help theoperator to quickly navigate the system and familiarize with all itsfeatures. On-line help is also available. A main screen for the suiteappears and displays, on its left side, several menus for a) selectingdifferent mathematical models, b) availability of experimental data, c)settings for graphical and tabular display of results, and d) buttonsfor starting a simulation experiment, stopping it, exiting the suitegoing back to the original suite selection screen, and a progress bar tomonitor the experiment evolution once started. On the right topside ofthe same screen there is a table where the results of the experiment aredisplayed and, on the right bottom side, a graph depicts the sameresults. The table can be set to be refreshed in real-time while thesimulation experiment is undergoing, or to be only updated at the end.The table can also be expanded, contracted, printed, and hidden. Thegraph can be refreshed either in real-time or at the end of theexperiment, expanded, contracted, printed, depicted in two or threedimensions, and the scales in the three axes can be selected as linearor logarithmic. The 3D graph can be rotated at will in real-time toexamine any particular structure by simply dragging the mouse. Bydouble-clicking on the graph, all graph features can be changed to suitoperator's needs and appeal. Different mathematical models can beselected from a menu list of single and multiple-scattering modelsemployed for the last three or four decades (ECAH, Mie, Waterman &Truell, Twersky, etc.) with the purpose of comparing their accuracy withthe last option in the menu which is the model for the current invention(METAMODEL™). The first model (Lamb-Epstein-Urick-Flammer) is awell-known combination of explicit approximations to the solution of thesingle-scattering wave equations. The rest of the models arefundamental. Several models can be selected as part of a singleexperiment by checking all the desired ones. Right-clicking on allsingle-scattering models but the first one opens a dialog window thatallows the operator to either select full-convergence or choose themaximum number of terms to be used in the series expansion of the fields(the range for the index n in the specification for the predictionengine). A recommended minimum for accuracy is also shown on the dialogscreen. This feature permits the researcher to understand theconvergence intricacies for different regimes, particle size, frequency,etc. Right clicking on all multiple-scattering models but the METAMODEL™one opens a dialog window that, besides the convergence feature, allowsthe operator to select the medium/particle signature, e.g. the Mie Modelif in the electromagnetics suite. Right clicking on the METAMODEL™model, a dialog screen appears that, when the control button is set tomanual, entitles the researcher to set, besides the medium/particlesignature and the convergence criterion, the following parameters forthe pair-correlation function: a) amplitude, b) oscillation frequency,and c) vanishing rate. For each one of these variables the user can seta minimum, a maximum, and a number of values in the range. Upon startingthe experiment, predictions for each one of those values will beconducted and results stored for post-processing and viewing.

[0087] Hiding the table unveils three menu frames. The one at the tophas two command buttons, one for specifying the nature of the phases ofthe composite, and the other for specifying the particle sizedistribution and concentration of the particulate phase. The two framesat the bottom permit the user to set the minimum, maximum, number ofvalues, and scale for the wave frequency (wavelength if in theelectromagnetic suite) as well as for the scattering angle. The legendon the command button at the left side of the top frame corresponds tothe nature of the phases (e.g. “Titania in Water”) and, when clicking onit a dialog window appears on the screen. Up to two different types ofparticles in one type of host medium can be specified from a database ofmaterials. Physical properties for the selected materials are displayedin a folder-like tab-strip control. New types of materials can be easilycreated with their physical properties—in a selected system of physicalunits—typed by the researcher and stored into a database of materials.Minimum, maximum, scale and number of values in the range to be use inthe simulation experiment can also be selected. Clicking on the SizeDistribution/Concentration command button pops up a dialog screen forsetting the size distribution and concentration of the particulatephase. Minimum, maximum, scale and number of values for the particleconcentration can be set by the operator. Regarding the particle sizedistribution, two populations can be specified (each made up of the sameor different particle material) and mathematically combined. Eachpopulation distribution can be modeled using a lognormal distribution.Minimum, maximum, scale, and number of values for the geometric mean anddeviation parameters for each population can be selected. As anergonomic feature to familiarize the user with this modeling of sizedistributions, three real images depicting special cases of a monosizedistribution, a wide unimodal distribution, and a bimodal distributionare shown at the top this screen and, when clicking on them thedistribution parameters are automatically set and displayed. Otherdistribution models besides the lognormal for each population canobviously be added to the interface if necessary to represent actualparticle size distributions for special cases. Clicking on the availableoption button in the Experimental Data frame triggers the appearance ofa dialog screen where: a) existing experimental data can be retrieved indifferent formats (binary, ASCII text files, databases, etc); b)Existing experimental data can be typed into the computer and storedinto its data base as a text file. Having completely defined asimulation experiment by specifying all the attributes of theconstituent phases, their ranges and discrete values, and the operatingvariables, the operator can start the experiment by clicking on theexecute button at the bottom left corner of the main screen. Theexperiment can be very complex and demanding on computer resourcesdepending upon the number of attributes that are to be varied and thenumber of points for each of them. Basic combinatorial calculations aredone by the system and displayed with the progress bar so that theoperator knows the status of the simulation experiment. Several thousandpredictions can be easily performed in a typical experiment. Before,during, or at the end of an experiment, clicking on the Graph/TableSettings button triggers a menu to appear overlapping the MathematicalModels and Experimental Data menus. These settings decide what is itthat the researcher wants to look at; it does not decide the nature(attributes to be changed and their ranges) of the simulationexperiment. At its top the experimenter can select one of fourattribute-axes (X, Y, Z, H). Clicking on one of them pops up a menu ofphysical attributes to be chosen for the particular axis. The Z-axis isused to display the attributes of the composite; the X, Y, and H axesare used to depict different attributes for the particulate and hostphases as well as the operating conditions. In the preferred embodimentbeing described, there are fifty-four different physical attributes tobe selected from for each of these three axes. If the graph is set astwo-dimensional, the plot will depict Z as a function of X parameterizedwith Y; if three-dimensional the plot will show Z as a function of X andY. The H-axis (Hyper-attribute) is a fourth-dimension attribute whicheffect on the simulation results can be graphically observed by varyingits value using a sliding rule control after having clicked on theH-Attribute button, selected the desired attribute, and dragging on thesliding control. Table and graph will dynamically update as the operatorslides the ruler through the range for the particular attributeselected. Clicking on the Z-axis button, displays a menu of physicalattributes for the composite to be shown on that axis. Among them are:attenuation, velocity, intensity, scattering cross-section, physicalproperties of the composite, extinction efficiency, etc. All thesealternative ways of viewing the relationship between differentattributes can be pursued in real-time while a complex simulationexperiment is undergoing.

SCOPE OF INVENTION

[0088] While my above description contains much specificity, that shouldnot be construed as limiting the scope of the invention, but rather asan exemplification of one preferred embodiment thereof Many othervariations are possible. Accordingly, the scope of the invention shouldbe determined not by the embodiment illustrated, but by the appendedclaims and their legal equivalents.

REFERENCES U.S. PATENT DOCUMENTS

[0089] [P1] U.S. Pat. No. 5,121,629, Jun. 16, 1992, Alba, F.

[0090] [P2] U.S. Pat. No. 5,619,324, Apr. 08, 1997, Harvill, T. L. andHolve, D. J.

[0091] [P3] U.S. Pat. No. 5,818,583, Oct. 06, 1998, E. Sevick-Muraca etal.

[0092] [P4] U.S. Pat. No. 6,109,098, Aug. 29, 2000, Dukhin, A. andGoetz, P.

[0093] [P5] U.S. Pat. No. 6,119,510, Sep. 19, 2000, Carasso, M. et al

COMMERCIAL ADVERTISEMENTS

[0094] [C1] MATLAB, MathWorks, Inc., Natick, Mass.

[0095] [C2] Maple V®, Waterloo Maple Software (800) 633-3345.

[0096] [C3] SCIENTIFIC WORKPLACE™, TCI Software Research, (800)622-3345.

[0097] [C4] MathCad 2001, MathSoft, (800) 628-4223.

[0098] [C5] Mathematica, Wolfram Research, Inc., (800) 622-3345

[0099] [C6] IMSL 32-bit Numerical Libraries and Presentation Graphicsfor Fortran and C, Visual Numerics, (713) 954-6785.

[0100] [C7] FEMLAB®, COMSOL, Inc., Burlington, Mass., (781) 273-3322.

[0101] [C8] Max-i, A Visual Electromagnetics Platform, Swiss FederalInstitute of Technology, Zurich Switzwerland.

[0102] [C9] Real Time, A Two-Dimensional Electromagnetic FieldSimulator, CRC Press, (800) 272-7737.

SCIENTIFIC PUBLICATIONS

[0103] [1] J. A. Stratton, “Electromagnetic Theory”, McGraw-Hill BookCompany, Inc. 1941.

[0104] [2] L. L. Foldy, “The Multiple Scattering of Waves—I. GeneralTheory of Isotropic Scattering by Randomly Distributed Scatterers”,Physical Review, Volume 67, Numbers 3 and 4, February 1 and 15, 1945.

[0105] [3] M. Lax, “Multiple Scattering of Waves”, Reviews of ModernPhysics, Volume 23, Number 4, October 1951.

[0106] [4] M. Lax, “Multiple Scattering of Waves. II. The EffectiveField in Dense Systems”, Physical Review, Volume 85, Number 4, Feb. 15,1952.

[0107] [5] P. S Epstein and R. R. Carhart, “The Absorption of Sound inSuspensions and Emulsions, I. Water Fog in Air”, J. Acoust. Soc. Am., 25[3] 553-565 (1953).

[0108] [6] A. R. Edmonds, “Angular Momentum in Quantum Mechanics”,Princeton University Press, 1957.

[0109] [7] N. G. Einspruch and R. Truell, “Scattering of a PlaneLongitudinal Wave by a Spherical Fluid Obstacle in an Elastic Medium”,J. Acoust. Soc. Am., 32 [2], 214-220, February 1960.

[0110] [8] N. G. Einspruch et al, “Scattering of a Plane Transverse Waveby a Sherical Obstacle in an Elastic Medium”, J. Appl. Phys., 31 [5],806-817, May 1960.

[0111] [9] P. C. Waterman & R. Truell, “Multiple Scattering of Waves”,Journal of Mathematical Physics, Volume 2, Number 4, July-August, 1961.

[0112] [10] V. Twersky, “On a General Class of Scattering Problems”, J.Mathematical Physics, 3 [11], July-August 1962.

[0113] [12] V. Twersky, “Multiple Scattering by Arbitrary Configurationsin Three Dimensions”, J. Mathematical Physics, 3 [1], January-February,1962.

[0114] [13] J. G. Fikioris & P. C. Waterman, “Multiple Scattering ofWaves. “Hole Corrections” in the Scalar Case”, Journal of MathematicalPhysics, Volume 5, Number 10, October 1964.

[0115] [14] J. C. F. Chow, “Attenuation of Acoustic Waves in DiluteEmulsions and Suspensions”, J. Acoust. Soc. Am., 36 [12] 2395-2401(1964).

[0116] [15] N. C. Mathur & K. C. Yeh, “Multiple Scattering ofElectromagnetic Waves by Random Scatterers of Finite Size”, Journal ofMathematical Physics, Volume 5, Number 11, November 1964.

[0117] [16] V. Twersky, “Acoustic Bulk Parameters of Random VolumeDistributions of Small Scatterers”, J. Acoust. Soc. Am., 36 [7], July1964.

[0118] [17] P. Lloyd and M. V. Berry, “Wave propagation through anassembly of spheres IV. Relations between different multiple scatteringtheories”, Proc. Phys. Soc., 1967, Vol. 91, [18] R. O. Watts and D.Henderson, “Pair Distribution Function for Fluid Hard Spheres”,Molecular Physics, 16 (3), 217-233, 1969.

[0119] [19] J. R. Allegra and S. A. Hawley, “Attenuation of Sound inSuspensions and Emulsions: Theory and Experiments”, J. Acoust. Soc. Am.,51, 1545-1564 (1972).

[0120] [20] A. K. Mal and S. K. Bose, “Dynamic Elastic Moduli of aSuspension of Imperfectly Bonded Spheres”, Proc. Camb. Phil. Soc. 76,587-600, 1974.

[0121] [21 ] V. Twersky, “Coherent Scalar Field in Pair-CorrelatedRandom Distributions of Aligned Scatterers”, J. Mathematical Physics, 18[12], December 1977.

[0122] [22] A. Ishimaru, “Wave Propagation and Scattering in RandomMedia”, Academic Press, 1978.

[0123] [23] V. Twerky, “Coherent Electromagnetic Waves inPair-Correlated Random Distributions of Aligned Scatterers”, J. Math.Phys. 19 (1), January 1978.

[0124] [24] A. J. Devaney, “Multiple Scattering Theory for Discrete,Elastic, Random Media”, J. Math. Phys. 21 9110, November 1980.

[0125] [25] V. N. Bringi et al, “Bulk Propagation Characteristics ofDiscrete Random Media”, Multiple Scattering and Waves in Random Media.P. L. Chow, W. E. Kohler, G. C. Papanicolau (eds.). North-HollandPublishing Company, 1981.

[0126] [26] L. Tsang and J. A. Kong, “Multiple Scattering of AcousticWaves by Random Distributions of Discrete Scatterers with the use ofQuasicrystalline-Coherent Potential Approximation”, Journal of AppliedPhysics (J.Appl.Phys.), 52(9), September 1981, pp. 5448-5458.

[0127] [27] van de Hulst, H. C., “Light Scattering by Small Particles”,Dover Publications, Inc., New York, 1981.

[0128] [28] L. Tsang et al, “Multiple Scattering of Acoustic Waves byRandom Distribution of Discrete Scatterers with the Quasicrystalline andPercus-Yevick Approximation”, Journal of the Acoustical Society ofAmerica (J. Acoust. Soc. Am.), 71(3), March 1982, pp 552-558.

[0129] [29] V. N. Bringi et al, “The Effects of Pair CorrelationFunction on Coherent Wave Attenuation in Discrete Random Media”, IEEETransactions on Antennas and Propagation, Vol. AP-30, No. 4, July 1982,pp. 805-808.

[0130] [30] L. Tsang and J. A. Kong, “Effective Propagation Constantsfor Coherent Electromagnetic Wave Propagation in Media Embedded withDielectric Scatterers”, J. Appl. Phys. 53(11), November 1982, pp.7162-7173.

[0131] [31] V. N. Bringi et al, “Coherent Wave Attenuation by a RandomDistribution of Particles”, Radio Science, Volume 17, Number 5, pages946-952, September-October 1982.

[0132] [32] A. Ishimaru and Y. Kuga, “Attenuation constant of a coherentfield in a dense distribution of particles”, Journal of the OpticalSociety of America (J. Opt. Soc. Am.), Vol. 72, No. 10, October 1982.

[0133] [33] V. K. Varadan et al, “Multiple Scattering Theory for Wavesin Discrete Random Media and Comparison with Experiments”, RadioScience, Volume 18, Number 3, pages 321-327, May-June 1983.

[0134] [34] V. K. Varadan et al, “Coherent Attenuation of Acoustic Wavesby Pair-Correlated Random Distributions of Scatterers with Uniform andGaussian Size Distributions”, J. Acoust. Soc. Am. 73(6), pp. 1941-1947,June 1983.

[0135] [35] Y. Ma et al, “Application of Twersky's Multiple ScatteringFormalism to a Dense Suspension of Elastic Particles in Water”, J.Acoust. Soc. Am. 75(2), February 1984, pp. 335-339.

[0136] [36] V. K. Varadan et al, “A Multiple Scattering Theory forElastic Wave Propagation in Discrete Random Media”, J. Acoust. Soc. Am.77(2), February 1985, pp. 375,385.

[0137] [37] V. V. Varadan et al, “Multiple Scattering Theory forAcoustic, Electromagnetic, and Elastic Waves in Discrete Random Media”,Multiple Scattering of Waves in Random Media and Random Rough Surfaces.The Pennsylvania State University, 1985.

[0138] [38] L. M. Schwartz, “Classical and Quantum Mechanical WavePropagation in Close Packed Disordered Suspensions”, Multiple Scatteringof Waves in Random Media and Random Rough Surfaces, The PennsylvaniaState University, 1985.

[0139] [39] W. A. Steele, “Particle Distributions at High Density”,Multiple Scattering of Waves in Random Media and Random Rough Surfaces,The Pennsylvania State University, 1985.

[0140] [40] E. D. Hirleman, “Modeling of Multiple Scattering Effects inFraunhofer Diffraction Particle Size Analysis”, ParticleCharacterization 5, 57-65 (1988).

[0141] [41] E. D. Hirleman, “A General Solution to the InverseNear-Forward Scattering Particle Sizing Problem in Multiple ScatteringEnvironments: Theory”, Proceedings of the 2^(nd) International Congresson Optical Particle Sizing, Mar. 5-8 1990, pp. 159-168.

[0142] [42] Y. Ma et al, “Comments on ultrasonic propagation insuspensions”, Journal of the Acoustical Society of America (J. Acoust.Soc. Am.) 87(6), June 1990.

[0143] [43] D. D. Phanord and N. E. Berger, “Multiple Scattering ofElastic Waves by a Distribution of Identical Spheres”, J. Math. Phys. 31(11), November 1990.

[0144] [44] J. Y. Kim et al, “Dispersive Wave Propagation in aViscoelastic Matrix Reinforced by Elastic Fibers”, J. Acoust. Soc. Am.,95 (3), March 1994.

[0145] [45] de Daran, F. et al, “Modeling of Electromagnetic WavesScattered by a system of spherical particles”, IEEE Transactions onMagnetics, Vol. 31, No. 3, May 1995.

[0146] [46] E. Sevick-Muraca et al, “Photon-Migration Measurement ofLatex Size Distribution in Concentrated Suspensions”, AIChE Journal,March 1997, Vol. 43, No. 3.

[0147] [47] C. Verdier and M. Piau, “Acoustic Wave Propagation inTwo-phase viscoelastic fluids: The Case of Polymer emulsions”, J.Acoust. Soc. Am. 101 [4], 1868-1876, April 1997.

[0148] [48] F. Alba et al, “Ultrasound Spectroscopy: A Sound Approach toSizing of Concentrated Particulates”, Handbook on Ultrasonic andDielectric Characterization Techniques for Suspended Particulates, TheAmerican Ceramic Society, 1998).

I claim:
 1. A computer-based method for allowing a human operator toconduct complex simulation, parametric studies, analysis and synthesisof real systems in which there exists interaction between physical wavesof arbitrary nature and a composite of particulate phases included in acontinuous phase, with arbitrary physical attributes, comprising thesteps of: predicting (21) the attributes of said composite from theknowledge of the attributes of the constituent phases, using afundamental mathematical model, to thereby obtain predictions that canbe contrasted with actual experimental values; inverting (22) saidmathematical model employing numerical optimization techniques whereby aselected subset of said attributes of said constituent phases andcomposite can be estimated from the knowledge of the rest of saidattributes of said constituent phases and composite; allowing said humanoperator to ergonomically interact (23) with said computer, whereby saidoperator can efficiently set up complex simulation experiments to assistin parametric studies, analysis, and synthesis of real systems in whichphysical waves of arbitrary nature interact with particulate composites.2. A method as set forth in claim 1 wherein said fundamentalmathematical model is based on the analytical wave scattering theory. 3.A method as set forth in claim 1 wherein said fundamental mathematicalmodel provides the means for successfully integrating viscous-inertial,thermal diffusion, surface tension, and viscoelasticity phenomena intothe fundamental statistical treatment of multiple-scattering of saidphysical waves interacting with said composite.
 4. A method as set forthin claim 1 wherein said predicting step further includes the use ofother mathematical models already existing in the prior-art whereby saidhuman operator can assess the differences in accuracy for each one ofsaid mathematical models when compared with actual experimental data. 5.A method as set forth in claim 1 wherein said composite has a particleconcentration in excess of 10% by volume.
 6. A method as set forth inclaim 1 wherein said composite comprises of solid inclusions in a solidhost medium.
 7. A method as set forth in claim I wherein said compositecomprises of solid particles suspended in a fluid host medium.
 8. Amethod as set forth in claim 1 wherein said composite comprises of fluidparticles suspended in a fluid host medium.
 9. A method as set forth inclaim 1 wherein said composite comprises of fluid particles suspended ina solid host medium.
 10. A software-implemented fundamental method ofpredicting the physical attributes of a highly concentrated composite ofparticulate phases with known attributes included in a continuous phasewith known attributes, in their interaction with physical waves ofarbitrary nature, comprising the steps of: characterizing statisticallythe field magnitudes in said composite by averaging said fieldmagnitudes throughout the ensemble of all possible spatial distributionsof particles and all possible distributions of particle attributes (24);characterizing deterministically the relation between the exciting andscattered fields around a particle with a common mathematical formalismregardless of the physical nature of said fields (25); characterizingstatistically the local spatial and attribute correlations betweenparticles (26); expressing the expected exciting field around a fixedparticle in terms of the incident field and the scattered fieldsproduced by the rest of the particles in the ensemble (27); averagingthroughout all possible attribute distributions of the said fixedparticle to thereby arrive to an infinite set of integral equationsrelating the expected fields around all particles (28); calculating saidintegrals, truncating said infinite set of equations, and specifying theoperating conditions to thereby arrive to a finite set of homogeneousequations in the field amplitudes (29); finding the complex root for thedeterminant of said homogeneous system of equations thereby to determinethe expected propagation constant for said composite (29); calculatingexpected values of said composite attributes from said expectedpropagation constant and the known attributes of said particulate phasesand said continuous phase (30).
 11. A method as set forth in claim 10wherein electromagnetic waves interact with suspensions of solidparticles in a fluid host medium with particle concentrations in excessof 1% by volume.
 12. A method as set forth in claim 10 whereinelectromagnetic waves interact with emulsions of fluid particles in afluid host medium with particle concentrations in excess of 1% byvolume.
 13. A method as set forth in claim 10 wherein acoustic wavesinteract with suspensions of solid particles in a fluid host medium withparticle concentrations in excess of 10% by volume.
 14. A method as setforth in claim 10 wherein acoustic waves interact with emulsions offluid particles in a fluid host medium with particle concentrations inexcess of 10% by volume.
 15. A method as set forth in claim 10 furtherincluding means for interfacing its software implementation withexisting acoustic and optical particle size measurement instruments withthe purpose of increasing the concentration range of said instruments.16. An apparatus for allowing a human operator to conduct complexsimulation, parametric studies, analysis and synthesis of real systemsin which there exists interaction between physical waves of arbitrarynature and a composite of particulate phases included in a continuousphase, with arbitrary physical attributes, comprising: a general-purposecomputer (19); data-carrier media (20) to store, distribute, and installsoftware and experimental data into said computer; predicting means(21), in combination with said computer and said carrier media, forpredicting the attributes of said composite from the knowledge of theattributes of said particulate phases and said continuous phase, using afundamental mathematical model, to thereby obtain predictions that canbe contrasted with actual experimental values; inverting means (22), incombination with said computer and carrier media, for inverting saidmathematical model using numerical optimization techniques, whereby aselected subset of said attributes of said particulate phases, saidcontinuous phase, and said composite can be estimated from the knowledgeof the rest of said attributes of said particulate phases, saidcontinuous phase, and said composite; interfacing means (23), incombination with said computer and said carrier media, for allowing saidhuman operator to interact ergonomically with said computer through agraphical interface, whereby said operator can efficiently set upcomplex simulation experiments to assist in parametric studies,analysis, and synthesis of real systems in which physical waves of anynature interact with particulate composites.
 17. An apparatus as setforth in claim 16 wherein the computer programming for said predictingand inverting means is implemented as a dynamically linked library(DLL).
 18. An apparatus as set forth in claim 16 wherein the computerprogramming for said predicting and inverting means is implementedthrough a Component Object Model (COM) server.
 19. An apparatus as setforth in claim 16, wherein said data-carrier media (20) comprises of atleast one compact disk or equivalent to store, distribute, and installsaid software and experimental data into said computer (19).
 20. Anapparatus as set forth in claim 16, wherein said data-carrier media (20)comprises storing, distributing, and installing said software andexperimental data into said computer (19) through the Internet.